寫(xiě)在最前面

本次的GATK實(shí)戰(zhàn)用到的raw data是水稻的宾尚，reference是v4版本泣崩，具體如下：
raw data：EV19-0_HQ_clean_R1.fq/EV19-0_HQ_clean_R2.fq
reference：IRGSP_v4.fasta
我在運(yùn)行GATK的時(shí)候并沒(méi)有每一步都用GATK的工具膀息，還涉及到另外幾個(gè)軟件，名稱和版本如下：
bwa-0.7.17
samtools-1.8
picard-tools-1.119
GenomeAnalysis TK-3.4.-0
將這些軟件安裝好并且加入到環(huán)境變量之后就可以進(jìn)入實(shí)戰(zhàn)部分啦。

以下是我的GATK實(shí)戰(zhàn)內(nèi)容：

1. 建立ref的index

bwa index IRGSP_v4.fasta

用bwa的index命令 建立reference序列的index芯肤，此步驟結(jié)束會(huì)產(chǎn)生5個(gè)文件，分別為IRGSP_v4.fasta為前綴的sa压鉴、pac崖咨、bwt、ann油吭、amb格式的文件

2. 建立ref的索引內(nèi)容

samtools faidx IRGSP_v4.fasta

用samtools的faidx命令 建立reference序列的索引內(nèi)容击蹲，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)以IRGSP_v4.fasta為前綴的fai文件

3. 建立ref的dict內(nèi)容

java -Xmx45g -jar CreateSequenceDictionary.jar R=IRGSP_v4.fasta O=IRGSP_v4.dict

用picard工具包 中的CreateSequenceDictionary建立reference序列的dict，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)以IRGSP_v4.fasta為前綴的dict文件

注：前三步結(jié)束一共會(huì)產(chǎn)生以ref名為前綴的7個(gè)文件婉宰，在此基礎(chǔ)上才能進(jìn)行后續(xù)命令程序歌豺，若在之前已經(jīng)進(jìn)行過(guò)以此ref為基礎(chǔ)的相關(guān)內(nèi)容，則可以跳過(guò)前三步

4. 對(duì)raw data文件進(jìn)行過(guò)濾心包，加header并生成sam文件

bwa mem -M -t 20 -R "@RG\tID:EV\tPL:ILLUMINA\tLB:EV\tSM:EV"  IRGSP_v4.fasta EV19-0_HQ_clean_R1.fq EV19-0_HQ_clean_R2.fq -o EV.sam

用bwa的mem命令 對(duì)fq文件進(jìn)行比對(duì)過(guò)濾类咧，并加header生成sam文件

參數(shù)解釋：
-M：用來(lái)處理同一個(gè) reads比對(duì)到參考基因組上不同位置的情況
-t：線程數(shù)
-R：設(shè)置 reads的 header
-o：輸出文件

注1：關(guān)于加header的內(nèi)容，一般格式為“@RG ID:樣品ID號(hào) PL：測(cè)序平臺(tái)（一般為ILLUMINA） LB：樣品的文庫(kù)名（一般為樣品名） SM：樣品名稱”其中所有的空格都是\t（制表符）
注2：對(duì)于pair-end文件，命令如上所示轮听，在ref后輸入 read1.fq read2.fq 骗露；對(duì)于sngle-end文件，則只需要在ref后輸入read.fq

5. 將sam文件轉(zhuǎn)為bam文件

samtools view -bS EV.sam -o EV.bam

用samtools的view命令 將sam文件轉(zhuǎn)為bam文件血巍，（bam文件是二進(jìn)制萧锉，運(yùn)算速度會(huì)快），此步驟結(jié)束會(huì)產(chǎn)生一個(gè)bam文件

參數(shù)解釋：
-b：輸出bam文件（因?yàn)槟J(rèn)輸出是sam）
-S：輸入sam文件（因?yàn)槟J(rèn)輸入是bam）

6. 對(duì)bam文件進(jìn)行排序

java -Xmx45g -jar SortSam.jar INPUT=EV.bam OUTPUT=EV_sort.bam SORT_ORDER=coordinate

用picard工具包 中的SortSam對(duì)上一步得到的bam文件中同一染色體對(duì)應(yīng)的條目按照坐標(biāo)順序從小到大進(jìn)行排序述寡，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)排序過(guò)的bam文件柿隙，可命名為EV_sort.bam

7. 對(duì)排序后的bam文件進(jìn)行去重

java -Xmx100g -jar MarkDuplicates.jar INPUT=EV_sort.bam OUTPUT=EV_dedup.bam METRICS_FILE=EV_dedup.metrics ASSUME_SORTED=true VALIDATION_STRINGENCY=SILENT CREATE_INDEX=true

用picard工具包 中的MarkDuplicates對(duì)排序后的bam文件進(jìn)行去重（由于PCR擴(kuò)增底物的過(guò)程中使樣品出現(xiàn)重復(fù)，故要去重鲫凶，此步使重復(fù)的reads帶上標(biāo)簽禀崖，便于后續(xù)處理），此步驟結(jié)束會(huì)產(chǎn)生一個(gè)metrics文件螟炫，一個(gè)bai文件波附，一個(gè)去重的bam文件，可命名為EV_dedup.bam
參數(shù)解釋：
ASSUME_SORTED：默認(rèn)為false昼钻，經(jīng)過(guò)sortsam后改為true
CREATE_INDEX：創(chuàng)建索引（true才會(huì)產(chǎn)生bai文件）

注1：metrics文件用于寫(xiě)入一些duplicates的統(tǒng)計(jì)信息
注2：CREATE_INDEX使得bai文件產(chǎn)生掸屡，若不在這一步創(chuàng)建，則下一步要先用samtools的index命令對(duì)EV_dedup.bam建立索引才行然评。

8. 建立indel區(qū)間文件

java -Xmx45g -jar GenomeAnalysisTK.jar -T RealignerTargetCreator -R IRGSP_v4.fasta -I EV_dedup.bam -o EV_dedup.intervals

用GATK工具包 中的RealignerTargetCreator先建立indel區(qū)間文件仅财，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)intervals文件
參數(shù)解釋：
-T：使用的工具名

9. 對(duì)以上所得位置進(jìn)行重新對(duì)比

java -Xmx45g -jar GenomeAnalysisTK.jar -T IndelRealigner -R IRGSP_v4.fasta -I EV_dedup.bam -targetIntervals EV_dedup.intervals -o EV_dedup_realign.bam

用GATK工具包 中的IndelRealigner對(duì)上一步建立的indel區(qū)域進(jìn)行reads重比對(duì)（indel會(huì)導(dǎo)致附近的錯(cuò)配，所以需要借由這一步降低indel附近的假陽(yáng)性）碗淌，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)重新比對(duì)過(guò)的bam文件盏求，可命名為EV_dedup_realign.bam

10. 進(jìn)行變異檢測(cè)

java -Xmx45g -jar GenomeAnalysisTK.jar -T UnifiedGenotyper -R IRGSP_v4.fasta -I EV_dedup_realign.bam -glm BOTH -stand_call_conf 30 -stand_emit_conf 10 -o LD.vcf

用GATK工具包 中的UnifiedGenotyper對(duì)上述重新比對(duì)過(guò)的bam文件進(jìn)行變異檢測(cè)，尋找variant亿眠，此步驟結(jié)束會(huì)產(chǎn)生一個(gè)vcf文件和一個(gè)idx文件

參數(shù)解釋：
-glm：選擇檢測(cè)的變異類型碎罚，“SNP、INDEL缕探、BOTH”皆可選魂莫，默認(rèn)為SNP还蹲，采用BOTH表示同時(shí)檢測(cè)兩者
-stand_call_conf：在變異檢測(cè)過(guò)程中用于區(qū)分低質(zhì)量變異位點(diǎn)和高質(zhì)量變異位點(diǎn)的閾值爹耗，高于其會(huì)被視為高質(zhì)量，低于其會(huì)標(biāo)注LowQual
-stand_emit_conf：在變異檢測(cè)過(guò)程中谜喊，所容許的最小質(zhì)量值潭兽。只有大于等于這個(gè)設(shè)定值的變異位點(diǎn)會(huì)被輸出到結(jié)果中。

寫(xiě)在最后面

1. 運(yùn)行命令時(shí)斗遏，要注意使用的軟件山卦、輸入文件以及輸出文件的路徑，最好采用絕對(duì)路徑的命令輸入模式
2. 由于過(guò)程中會(huì)產(chǎn)生很多文件诵次，所以輸出文件要好好命名账蓉。每個(gè)文件是從哪一步得到的枚碗，最好標(biāo)上后綴，方便后續(xù)調(diào)用
3. 我這里的GATK版本還是較舊的3.4版本铸本，現(xiàn)在最新的應(yīng)該是4.0版本肮雨。在4.0版本中去掉了RealignerTargetCreator、IndelRealigner箱玷、UnifiedGenotyper這幾個(gè)工具怨规，全部用HaplotypeCaller可以解決，具體的參數(shù)在運(yùn)行程序的時(shí)候help一下就可以看到了锡足，所以8 波丰、9、10三步驟可以用HaplotypeCaller代替
4. 以上的GATK實(shí)戰(zhàn)步驟也是我在摸索中運(yùn)行的舶得，后續(xù)的數(shù)據(jù)分析還沒(méi)有進(jìn)行過(guò)掰烟，其中可能存在一定問(wèn)題，歡迎提出指正沐批，互相交流媚赖。