高橋鉄美｜フリーソフトを使った QTL 解析

ホーム > フリーソフトを使った QTL 解析

2020年12月24日作成
2021年1月2日修正
2021年6月17日R/qtl2の説明に共変量について追加

フリーソフトを使った QTL 解析

F0（祖父、祖母）とF2世代からddRAD-seqでデータを取得し、QTL解析をしました。備忘録としてまとめます。使用するソフトは、すべてフリーです。当ホームページによって生じた不利益については、一切の責任を負いかねますのでご了承ください。

1) Stacks: denovoアセブル
2) Lep-MAP3: 連鎖群（LGs）作成
3) MapChart: LGsの可視化
4) r/qtl2: 量的形質のマッピング

参考にしたページ
https://catchenlab.life.illinois.edu/stacks/manual/
https://avikarn.com/2019-04-17-Genetic-Mapping-in-Lep-MAP3/
https://kbroman.org/qtl2/assets/vignettes/user_guide.html

1) Stacks: denovoアセブル
※今回は、データがすでにフィルタリングされた状態で納品されたので、データを個体ごとに振り分けてフィルタリングするprocess_radtagsは実行しませんでした。
※各パイプラインの詳しい説明は、「Stacks の簡単な使い方」も参照してください。

1-1) ustacks: まず塩基配列が全く同じリードを束ね（stacks）、次にstacks同士を比較してSNPsを検出する。全ての個体（F0とF2の両方）で、１個体ずつ下のコマンドを実行する。

$ /Applications/stacks-2.54/ustacks -f ./rawdata/parent001.fastq.gz -o ./stacks/ -i 1 -m 3 -M 2 -p 3

-f: 解析するファイルを指定する（ここではペアエンドのうち、最初のファイルだけ解析、gz圧縮ファイルでも良い）。
-o: 出力ファイルを保存するフォルダ（ディレクトリ）を指定する。
-i: サンプルID。個体ごとに異なる数字を指定する。
-m: stacksを作る最低のカバレージ。カバレージがこれより小さいリードは、解析から排除される。
-M: ファイル内で同じオーソログとする最低距離（SNP数）。
-p: スレッド数

出力ファイル（個体ごとに３つのファイルが作られる）

parent001.alleles.tsv.gz
parent001.snps.tsv.gz
parent001.tags.tsv.gz
.
.
.

1-2) cstacks: F0個体間のSNPsを検出する。

F0の2個体（祖父、祖母）のみのpopulation mapを作成し、テキスト形式で保存する（popmap_parent.txt）。

parent001[タブ]parent[改行]
parent002[タブ]parent[改行]

※左のparent001などは、データファイル（parent001.fastq.gzなど）から拡張子（.fastq.gz）を省いた名前。ustacksの出力ファイルから拡張子を省いた名前とも一致する。
※右はparentとする。
※エクセルで作成してタブ区切りテキストで保存すると、なぜか上手くいかなかった。

下のコマンドを実行する。

$ /Applications/stacks-2.54/cstacks -P ./stacks/ -M ./popmap_parent.txt -n 2 -p 2

-P: 解析に使用するtsvファイル（ustacksの出力）が入っているフォルダを指定する。
-M: F0のみのpopulation mapファイルを指定する。
-n: 祖父と祖母の間で同じオーソログとする最低距離（SNP数）。
-p: スレッド数。

出力ファイル

catalog.alleles.tsv.gz
catalog.snps.tsv.gz
catalog.tags.tsv.gz

1-3) sstacks: 各個体のSNPsをジェノタイプする。

全個体を含んだpopulation mapを作成し、テキスト形式で保存する（popmap_all.txt）。

parent001[タブ]parent[改行]
parent002[タブ]parent[改行]
progeny003[タブ]progeny[改行]
progeny004[タブ]progeny[改行]
progeny005[タブ]progeny[改行]
.
.
.

下のコマンドを実行する。

$ /Applications/stacks-2.54/sstacks -P ./stacks/ -M ./popmap_all.txt -p 2

-P: 解析に使用するtsvファイル（ustacksとcstacksの出力）が入っているフォルダを指定する。
-M: 全個体を含んだpopulation mapを指定する。
-p: スレッド数。

出力ファイル（個体ごとに１つのファイルが作られる）

parent001.matches.tsv.gz
.
.
.

1-4) tsv2bam: ペアエンドのもう片方のリードをくっつけて、BAMファイルを作る。下のコマンドを実行する。

$ /Applications/stacks-2.54/tsv2bam -P ./stacks/ -M ./popmap_all.txt -R ./rawdata/ -t 2

-P: ustacks, cstacks, sstacksの出力ファイルの入ったフォルダを指定する。
-M: 全個体の入ったpopulation mapを指定する。
-R: ペアエンドの後ろの方（まだ解析していない方）のファイルが入ったフォルダを指定する。ペアエンドの前の方の名前が「hoge.fastq.gz」だった場合、後ろの方の名前はそれに「.2」を加えて「hoge.2.fastq.gz」にしておく。
-t: スレッド数。

出力ファイル例

parent001.matches.bam
.
.
.
tsv2bam.log

1-5) gstacks: 全体を通してSNPの特定とジェノタイプをする。下のコマンドを実行する。

$ /Applications/stacks-2.54/gstacks -P ./stacks/ -M ./popmap_all.txt -t 2

-P: ustacks, cstacks, sstacks, tsv2bamの出力ファイルの入ったフォルダを指定する。
-M: 全個体を含んだpopulation mapを指定する。
-t: スレッド数。

出力ファイル

catalog.calls
catalog.fa.gz
gstacks.log
gstacks.log.distribs

1-6) populations: vcfファイルを作る。

$ /Applications/stacks-2.54/populations -P ./stacks/ --vcf -t 2 -r 0.75 --write-single-snp

-P: 入出力ファイルのフォルダを指定する。
--vcf: SNPsをvcf形式で出力する。
-t: スレッド数。
-r: 必要最低カバレージ。ジェノタイピングできた個体の割合がこの値より高いローカスのみを出力する（低いローカスは削除される）。
--write-single-snp: ローカス内に複数のSNPsがある場合、最初のSNPのみを出力する。

出力ファイル

populations.haplotypes.tsv
populations.hapstats.tsv
populations.log
populations.log.distribs
populations.snps.vcf（SNPデータの入った出力ファイル）
populations.sumstats_summary.tsv
populations.sumstats.tsv

トップに戻る

2) Lep-MAP3: 連鎖群（LGs）作成
※javaを使用します。
※今回は、populations.snps.vcfをもとに、F0がaaxbbのローカスのみのファイル（aaxbb.snps.vcf）を作成して解析しました。

2-1) Pedigree fileをエクセルで作り、タブ区切りテキスト（pedigree.txt）で保存する。

CHR	POS	family_name	family_name	family_name	family_name	.....
CHR	POS	parent001	parent002	progeny003	progeny004	.....
CHR	POS	0	0	parent001	parent001	.....
CHR	POS	0	0	parent002	parent002	.....
CHR	POS	1	2	0	0	.....
CHR	POS	0	0	0	0	.....

１行目：上の通りにする。
２行目：parent001やprogeny003は、vcf入力ファイルで使っている個体名と一致させる。
３行目：全てのF2個体（progeny）について、祖父（F0オス）を指定。
４行目：全てのF2個体（progeny）について、祖母（F0メス）を指定。
５行目：性別を指定。祖父は1、祖母は2とする。F2個体は0でいいみたい。
６行目：意味はわからないが、全個体0とするみたい。

2-2) ParentCall2: 下のコマンドを実行する。

$ java -cp /Applications/Lep-MAP3/binary+code/bin ParentCall2 data = pedigree.txt vcfFile = aaxbb.snps.vcf > p.call

data: 上で作成したPedigree fileを指定する。
vcfFile: vcf入力ファイルを指定する。

2-3) Filtering2: 下のコマンドを実行する。

$ java -cp /Applications/Lep-MAP3/binary+code/bin Filtering2 data = p.call removeNonInformative = 1 dataTolerance = 0.001 > p_fil.call

data: 上の解析の出力ファイルを指定する。
removeNonInformative: 1を指定すると、情報のないローカスは削除される。
dataTolerance: 予想される遺伝子型分離比（この場合1:2:1）からのずれが許容される下限を指定する。

2-4) SeparateChromosomes2: ローカスをLGごとに分ける。下のコマンドを実行する。

$ java -cp /Applications/Lep-MAP3/binary+code/bin SeparateChromosomes2 data = p_fil.call lodLimit = 10 theta = 0.05 > map.txt

data: 上の解析の出力ファイルを指定する。
lodLimit: 同じLGにまとめる際に許容されるLODスコア。
theta: 組み替え率。

下のコマンドで、LGごとのローカス数を確認できる。

$ sort map.txt | uniq -c | sort -n
1 #java SeparateChromosomes2 data=p_fil.call lodLimit=10 theta=0.05
3 24
4 23
9 22
11 21
12 20
13 19
16 18
19 17
20 16
21 14
21 15
22 12
22 13
24 10
24 11
24 9
26 8
27 7
28 6
29 4
29 5
35 3
39 2
48 1

左列：ローカス数（map.txt内での行数）。
右列：LG番号。

2-5) OrderMarkers2: LGごとに、マーカーの位置（端からのcM）を計算する。下のコマンドを実行する。

$ java -cp /Applications/Lep-MAP3/binary+code/bin OrderMarkers2 data = p_fil.call map = map.txt sexAveraged = 1 > order01.txt

※計算ごとに異なるrandom seedが与えられ、結果が異なる。複数回（例えば20回）計算を繰り返し、LGごとに尤度が最大の結果をコピペして、一つのファイル（order.txt）にまとめる。

2-6) map2genotypes.awk: Lep-MAP3形式の出力ファイルを、一般的な出力ファイルに変換する。Avi Karnさんのホームページからmap2genotypes.awkをダウンロードし、上の解析の出力ファイルと同じフォルダに保存して、下のコマンドを実行する。

$ awk -v fullData=1 -f map2genotypes.awk order.txt > genotypes.txt

※「=（イコール）」の前後にスペースを入れない。
※出力ファイルgenotypes.txtのマーカー名が書き換えられているので、下の手順で元のマーカー名に書き換える。
1) エクセルでmap.txtを読み込み、1行目を削除する。
2) エクセルでp_fil.callを読み込み、A列のマーカー名（８行目以下）をコピーして、map.txtのB列にペーストする（A, B列でマーカー数が同じか確認）。
3) map.txtで、A列が「0」のマーカーを削除する（LGに振り分けられていないため）。
4) エクセルでgenotypes.txtを読み込み、A列の仮マーカー名（7行目以下）をコピーして、map.txtのC列にペーストする。
5) map.txtのD列に、上から1, 2, 3, ...と番号を振る。
6) map.txtでC, D列を選択し、C列でソート（最小から最大）する。
7) map.txtでA, B, C, D列を選択し、D列でソート（最小から最大）する（A列のLG番号が昇順で並んでいること、C列の仮マーカー名がgeotypes.txtのA列と同じ並びであることを確認）。
8) map.txtのB列のマーカー名をコピーし、genotypes.txtのA列7行目以下にペーストする（マーカー数が正しいか確認）。
9) map.txtのA列のLG番号をコピーし、genotypes.txtのB列7行目以下にペーストする。
10) genotypes.txtを保存する。

トップに戻る

3) MapChart: LGsの可視化
※MapChartのMac版はない。Windows版を使う。

3-1) エクセルを使って、genotypes.txtから下のファイルを作成し、タブ区切りテキスト（genotypes_MapChart.txt）で保存する。