GATA檢測(cè)變異超詳細(xì)教程

應(yīng)該是最詳細(xì)的生信分析教程，主要由新必應(yīng)（也就是搭載的chat GP T 4.0）協(xié)助完成，我只提供了分析代碼和做了最終調(diào)整排版。

ChatGPT對(duì)代碼進(jìn)行了逐行解釋和分節(jié)進(jìn)行總結(jié)說(shuō)明，同時(shí)對(duì)我的代碼進(jìn)行了糾錯(cuò)（確實(shí)是我沒(méi)有發(fā)現(xiàn)的小錯(cuò)誤，挺驚喜，以后讓AI進(jìn)行改腳本應(yīng)該效率挺高）

對(duì)于生信初學(xué)者絕對(duì)是很好的助攻，同樣對(duì)生信分析流程記錄和教程撰寫(xiě)來(lái)說(shuō)也是全新的體驗(yàn)，省去了很多對(duì)分析流程的調(diào)整和解釋說(shuō)明，使用起來(lái)非常方便，特別是對(duì)不熟悉的代碼的解讀以及對(duì)代碼進(jìn)行糾錯(cuò)

雖然只是AI帶來(lái)的編程學(xué)習(xí)和分析的一些便利，但總感覺(jué)這只是開(kāi)始，后面不知道會(huì)給我們的學(xué)習(xí)和生活帶來(lái)多少影響和顛覆

分析流程正文

1、定義變量和路徑

第一部分定義了一些變量，表示你后面要用到的路徑和文件


fasta=Malus_domestica_HFTH1_v1.0.a1.fa #定義參考基因組的文件名

TEMP=/data/temp #定義一個(gè)臨時(shí)目錄

PICARD=/data/software/picard.jar #定義picard工具的路徑

GATK=/data/Erick_Tong/software/gatk-4.1.0.0/ #定義gatk工具的路徑

WK=/data/02Data_analysis_project/01BSA_seq_analysis/HF_BSAseq #定義工作目錄

PK=/data/02Data_analysis_project/01BSA_seq_analysis #定義項(xiàng)目目錄

2、Trimmomatic清洗數(shù)據(jù)

  • 第二部分檢查是否有一個(gè)叫01data/clean的目錄，如果沒(méi)有就創(chuàng)建一個(gè)。這個(gè)目錄會(huì)存放去除接頭和低質(zhì)量堿基后的修剪和清洗過(guò)的序列。
  • 使用Trimmomatic工具來(lái)對(duì)sampleid.txt中列出的每個(gè)樣本的雙端序列進(jìn)行修剪和清洗。Trimmomatic需要一些參數(shù)來(lái)指定如何修剪序列，比如LEADING, TRAILING, SLIDINGWINDOW, MINLEN 和 TOPHRED33。
if [ ! -d ./01data/clean ] #如果沒(méi)有./01data/clean這個(gè)目錄

then mkdir -p ./01data/clean #就創(chuàng)建一個(gè)./01data/clean這個(gè)目錄

fi

for i in cat $WK/01data/sampleid.txt #對(duì)$WK/01data/sampleid.txt文件中列出的每個(gè)樣本名字做循環(huán)

#01rawdata

clean do trimmomatic PE -threads 12 \ #用trimmomatic工具對(duì)雙端序列進(jìn)行修剪和清洗，使用12個(gè)線程，并把結(jié)果輸出到下面四個(gè)文件中，分別是配對(duì)和不配對(duì)的序列。

$WK/01data/${i}_1.fq.gz $WK/01data/${i}_2.fq.gz \ #輸入原始雙端序列文件，${i}表示樣本名字。

$WK/01data/clean/${i}_1_paired_clean.fq.gz \ #輸出修剪和清洗后的第一條配對(duì)序列文件。

$WK/01data/clean/${i}_1_unpair_clean.fq.gz \ #輸出修剪和清洗后的第一條不配對(duì)序列文件。

$WK/01data/clean/${i}_2_paired_clean.fq.gz \ #輸出修剪和清洗后的第二條配對(duì)序列文件。

$WK/01data/clean/${i}_2_unpair_clean.fq.gz \ #輸出修剪和清洗后的第二條不配對(duì)序列文件。

LEADING:3 TRAILING:3 \ #指定去除每條序列開(kāi)頭和結(jié)尾質(zhì)量值低于3的堿基。

SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33 #指定使用一個(gè)長(zhǎng)度為4堿基、質(zhì)量值閾值為20的滑動(dòng)窗口來(lái)掃描每條序列，并去除窗口內(nèi)平均質(zhì)量值低于閾值的部分；指定去除長(zhǎng)度小于50堿基的序列；指定使用Phred+33編碼格式來(lái)表示質(zhì)量值。

done

3、clean data數(shù)據(jù)質(zhì)控

  • 檢查是否有一個(gè)叫01data/clean/fastqc的目錄，如果沒(méi)有就創(chuàng)建一個(gè)。這個(gè)目錄會(huì)存放修剪和清洗過(guò)的序列的質(zhì)量控制報(bào)告。

  • 使用一個(gè)叫FastQC 的工具來(lái)為每對(duì)修剪和清洗過(guò)的序列生成質(zhì)量控制報(bào)告。FastQC提供了一系列模塊化的分析方法來(lái)檢查序列質(zhì)量的各個(gè)方面，比如每個(gè)堿基位置上的質(zhì)量值、每條序列上的質(zhì)量值、每個(gè)堿基位置上堿基含量、每個(gè)堿基位置上GC含量、每條序列上GC含量、每個(gè)堿基位置上N含量、序列長(zhǎng)度分布、序列重復(fù)水平、過(guò)度表達(dá)序列等等。FastQC還生成HTML文件，包含圖表和表格來(lái)可視化質(zhì)量指標(biāo)。

  • 使用MultiQC工具來(lái)把FastQC報(bào)告匯總成一個(gè)總結(jié)報(bào)告，可以方便地瀏覽。

if [ ! -d ./01data/clean/fastqc ] #如果沒(méi)有./01clean/fastqc這個(gè)目錄

then mkdir -p ./01clean/fastqc #就創(chuàng)建一個(gè)./clean/fastqc這個(gè)目錄

fi

for i in cat $WK/01sampleid.txt #對(duì)$WK/sampleid.txt文件中列出的每個(gè)樣本名字做循環(huán)

do

fastqc $WK/data/clean/${i}_1_paired_clean.fq.gz -o $WK/data/clean/fastqc/ #用fastqc工具對(duì)修剪和清洗后的第一條配對(duì)序列生成質(zhì)量控制報(bào)告，并把結(jié)果輸出到$WK/data/clean/fastqc/

fastqc $WK/data/clean/${i}_2_paired_clean.fq.gz -o $WK/data/clean/fastqc/ #用fastqc工具對(duì)修剪和清洗后的第二條配對(duì)序列生成質(zhì)量控制報(bào)告，并把結(jié)果輸出到$WK/data/clean/fastqc/

done

multiqc $WK/data/clean/fastqc/ -o $WK/data/clean/fastqc/ #用multiqc工具對(duì)fastqc生成的報(bào)告進(jìn)行匯總，并把結(jié)果輸出到$WK/data/clean/fastqc/

4、比對(duì)

  • 進(jìn)行基因組比對(duì)的步驟，使用了bwa軟件包來(lái)實(shí)現(xiàn)。

  • 并進(jìn)行PCR重復(fù)序列的標(biāo)記和刪除的步驟，使用picard軟件包來(lái)實(shí)現(xiàn)。

  • 然后進(jìn)行比對(duì)結(jié)果的統(tǒng)計(jì)和評(píng)估的步驟，使用了samtools軟件包來(lái)實(shí)現(xiàn)。

if [ ! -d $WK/02mapping ] # 如果不存在目錄$WK/02mapping，就創(chuàng)建一個(gè)

then mkdir -p $WK/02mapping # 創(chuàng)建目錄$WK/02mapping

fi

bwa index -a bwtsw $REF/GDDH13_1-1_formatted2.fasta -p $REF/GDDH13 # 使用bwa index命令為參考基因組$REF/GDDH13_1-1_formatted2.fasta建立索引，索引文件保存在$REF/GDDH13中

for id in cat $WK/01data/sampleid.txt # 對(duì)于每一個(gè)樣本id，都執(zhí)行以下操作

do

bwa mem -t 12 -M -R “@RG\tID:$id\tSM:$id\tLB:WES\tPL:Illumina” \ # 使用bwa mem命令將樣本的雙端測(cè)序數(shù)據(jù)與參考基因組進(jìn)行比對(duì)，使用12個(gè)線程，標(biāo)記次優(yōu)比對(duì)位置，添加read group信息

$REF/GDDH13 \ # 指定參考基因組索引文件

$WK/01data/clean/${id}1_paired_clean.fq.gz \ # 指定第一條read文件

$WK/01data/clean/${id}2_paired_clean.fq.gz | samtools sort -@ 12 -o $WK/02mapping/bwa_mem${id}.bam # 指定第二條read文件，并將比對(duì)結(jié)果通過(guò)管道傳遞給samtools sort命令，使用12個(gè)線程將結(jié)果按照坐標(biāo)排序并輸出為bam格式文件，保存在$WK/02mapping/bwa_mem${id}.bam中java -Xmx64g -jar $PICARD MarkDuplicates \ # 使用java運(yùn)行picard的MarkDuplicates工具，分配64g內(nèi)存空間

I=$WK/02mapping/bwa_mem_${id}.bam \ # 指定輸入的bam文件為上一步的比對(duì)結(jié)果文件

O=$WK/02mapping/${id}_sort_rmdup.bam \ # 指定輸出的bam文件為去除重復(fù)后的文件，保存在$WK/02mapping/${id}_sort_rmdup.bam中CREATE_INDEX=true \ # 創(chuàng)建索引文件

REMOVE_DUPLICATES=true \ # 刪除重復(fù)序列而不只是標(biāo)記它們

M=$WK/02mapping/${id}_sort_rmdup_metrics.txt # 輸出重復(fù)序列的統(tǒng)計(jì)信息到文本文件中，保存在$WK/02mapping/${id}_sort_rmdup_metrics.txt中samtools flagstat $WK/02mapping/${id}_sort_rmdup.bam >
$WK/02mapping/${id}_deduped_stat.txt # 使用samtools flagstat命令對(duì)去除重復(fù)后的bam文件進(jìn)行統(tǒng)計(jì)分析，并將結(jié)果輸出到文本文件中，保存在$WK/02mapping/${id}_deduped_stat.txt中samtools coverage $WK/02mapping/${id}_sort_rmdup.bam > $WK/02mapping/${id}_deduped_coverage.txt # 使用samtools coverage命令計(jì)算去除重復(fù)后的bam文件的覆蓋度，并將結(jié)果輸出到文本文件中，保存在$WK/02mapping/${id}_deduped_coverage.txt中

done

5、檢測(cè)變異

  • 對(duì)每個(gè)樣本進(jìn)行HaplotypeCaller分析，并輸出GVCF格式的文件。這個(gè)步驟可以檢測(cè)出每個(gè)樣本中的變異位點(diǎn)和非變異位點(diǎn)，并將非變異位點(diǎn)分成塊以壓縮文件大小。

  • 用CombineGVCFs模塊將多個(gè)GVCF文件合并成一個(gè)大的GVCF文件，以便于后續(xù)進(jìn)行聯(lián)合分析。這個(gè)步驟可以將不同樣本中的變異信息整合到一起，并保留每個(gè)樣本的基因型似然值。

  • 用GenotypeGVCFs模塊對(duì)合并后的GVCF文件進(jìn)行聯(lián)合基因型分析，輸出最終的VCF文件。這個(gè)步驟可以根據(jù)所有樣本的數(shù)據(jù)計(jì)算出每個(gè)位點(diǎn)上每個(gè)樣本的最可能基因型，并給出質(zhì)量評(píng)估和過(guò)濾標(biāo)記。

#01檢測(cè)變異

if [ ! -d $WK/03vcf ] # 如果不存在目錄$WK/03vcf，就創(chuàng)建這個(gè)目錄

then mkdir -p $WK/03vcf # 創(chuàng)建目錄$WK/03vcf

fi # 結(jié)束if語(yǔ)句

for i in `cat $WK/01data/sampleid.txt` # 對(duì)$WK/01data/sampleid.txt文件中的每一行（即每個(gè)樣本ID）執(zhí)行循環(huán)

do # 開(kāi)始循環(huán)

gatk --java-options "-Xmx32g -Djava.io.tmpdir=$TEMP" HaplotypeCaller \ # 用gatk軟件包調(diào)用HaplotypeCaller模塊，設(shè)置java選項(xiàng)為-Xmx32g（最大內(nèi)存為為32G）和-Djava.io.tmpdir=$TEMP（臨時(shí)文件夾為$TEMP）
     -R $REF/GDDH13_1-1_formatted2.fasta \ # 設(shè)置參考基因組文件為$REF/GDDH13_1-1_formatted2.fasta
     -I $WK/02mapping/${i}_sort_rmdup.bam \ # 設(shè)置輸入文件為$WK/02mapping/${i}_sort_rmdup.bam，其中${i}是樣本ID
     -O $WK/03vcf/${i}_gvcf.gz \ # 設(shè)置輸出文件為$WK/03vcf/${i}_gvcf.gz，其中${i}是樣本ID，輸出格式為GVCF
     -ERC GVCF # 設(shè)置發(fā)射模式為GVCF，即只輸出變異位點(diǎn)和非變異位點(diǎn)的塊信息
 done # 結(jié)束循環(huán)

# 02合并多個(gè)GVCF文件
 ls $WK/03vcf/*_gvcf.gz > $WK/03vcf/gvcf_all.list # 列出所有以_gvcf.gz結(jié)尾的文件，并保存到$WK/03vcf/gvcf_all.list中
 /usr/bin/java -Xmx32G -jar $GATK/gatk-package-4.1.0.0-local.jar CombineGVCFs \ # 用/usr/bin/java軟件調(diào)用$GATK/gatk-package-4.1.0.0-local.jar包中的CombineGVCFs模塊，設(shè)置最大內(nèi)存為32G 
     -R $REF/$fasta \ # 設(shè)置參考基因組文件為$REF/$fasta 
    
 -V $WK/02HaplotypeCaller_gvcf/bwa_mem_A1_1_g_vcf.gz \ # 設(shè)置輸入文件之一為$WK/02HaplotypeCaller_gvcf/bwa_mem_A1_1_g_vcf.gz 
     -V $WK/02HaplotypeCaller_gvcf/bwa_mem_A2_gr_vcf.gz \ # 設(shè)置輸入文件之一為$WK/02HaplotypeCaller_gvcf/bwa_mem_A2_gr_vcf.gz 
     -V $PK/parents_genome/ss_genome/03vcf/ss_gvcf.gz \  # 設(shè)置輸入文件之一為$PK/parents_genome/ss_genome/03vcf/ss_gvcf.gz 
     -V $PK/parents_genome/hanfu1_genome/hanfu1_gvcf.gz \  # 設(shè)置輸入文件之一為$PK /parents_genome/hanfu1_genome/hanfu1_gvcf.gz 
     -O $WK /03Combine_vcf/A1A2_SSHF_Combine_G.vcf  # 設(shè)置輸出文件為$ WK /03Combine_vcf/A1A2_SSHF_Combine_G.vcf 

#03joint genotyping
 /usr/bin/java -Xmx32g -Djava.io.tmpdir=$TEMP -jar $GATK/gatk-package-4.1.0.0-local.jar GenotypeGVCFs \  # 用/usr/bin/java軟件調(diào)用$GATK/gatk-package-4.1.0.0-local.jar包中的GenotypeGVCFs模塊，設(shè)置最大內(nèi)存為32G和臨時(shí)文件夾為$TEMP 
        -R $REF/$fasta \  # 設(shè)置參考基因組
        -V $WK/03Combine_vcf/A1A2_SSHF_Combine_G.vcf \ # 設(shè)置輸入文件為$WK/03Combine_vcf/A1A2_SSHF_Combine_G.vcf，即合并后的GVCF文件
        -O $WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf  # 設(shè)置輸出文件為$WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf，即最終的VCF文件

6、 提取SNP（InDel）并進(jìn)行過(guò)濾

這一部分使用SelectVariants功能，從A1A2_SSHF_Genotype.vcf文件中提取出SNP類(lèi)型的變異，并輸出到01A1A2_SSHF_SNP_raw.vcf文件中。

#01 提取SNP

/usr/bin/java -Xmx32G -jar $GATK/gatk-package-4.1.0.0-local.jar SelectVariants \ 
-R $REF/$fasta \ 
-V $WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf \ 
--select-type SNP \ 
-O $WK/04Variants/01A1A2_SSHF_SNP_raw.vcf

#每行代碼的含義如下：

# `/usr/bin/java -Xmx32G`：使用java程序運(yùn)行GATK，并設(shè)置最大內(nèi)存為32GB。
# `-jar $GATK/gatk-package-4.1.0.0-local.jar`：指定運(yùn)行的GATK版本為4.1.0.0。
# `SelectVariants`：指定使用SelectVariants功能。
# `-R $REF/$fasta`：指定參考基因組序列文件，$REF和$fasta是變量，需要根據(jù)實(shí)際情況替換。
# `-V $WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf`：指定輸入的vcf文件，$WK是變量，需要根據(jù)實(shí)際情況替換。
# `--select-type SNP`：指定提取SNP類(lèi)型的變異。
# `-O $WK/04Variants/01A1A2_SSHF_SNP_raw.vcf`：指定輸出的vcf文件。
這一部分使用VariantFiltration功能，對(duì)01A1A2_SSHF_SNP_raw.vcf文件中的SNP進(jìn)行質(zhì)量過(guò)濾，并添加過(guò)濾標(biāo)簽到02A1A2_SSHF_SNP_raw_filter.vcf文件中。

#02 過(guò)濾SNP

/usr/bin/java -Xmx32G -jar $GATK/gatk-package-4.1.0.0-local.jar VariantFiltration \ 
-R $REF/$fasta \ 
-V $WK/04Variants/01A1A2_SSHF_SNP_raw.vcf \ 
--cluster-window-size 10 --cluster-size 3 \ 
--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || DP < 8.0 " \ 
--filter-name 'SNP_filter' \ 
--missing-values-evaluate-as-failing \ 
-O $WK/04Variants/02A1A2_SSHF_SNP_raw_filter.vcf

每行代碼的含義如下：

# `/usr/bin/java -Xmx32G`：同上。
# `-jar $GATK/gatk-package-4.1.0.0-local.jar`：同上。
# `VariantFiltration`：指定使用VariantFiltration功能。
# `-R $REF/$fasta`：同上。
# `-V $WK/04Variants/01A1A2_SSHF_SNP_raw.vcf`：指定輸入的vcf文件，與上一步輸出相同。
# `--cluster-window-size 10 --cluster-size 3`：指定過(guò)濾掉窗口大小為10bp內(nèi)有3個(gè)或以上變異位點(diǎn)的記錄。
# `--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || DP < 8.0 ：指定過(guò)濾表達(dá)式，根據(jù)一些質(zhì)量指標(biāo)對(duì)變異進(jìn)行過(guò)濾。這些指標(biāo)包括：
#QD：每個(gè)變異位點(diǎn)的質(zhì)量值除以深度。
#MQ：比對(duì)到參考基因組的平均質(zhì)量值。
#FS：Fisher’s exact test檢測(cè)是否有分層偏差。
#SOR：StrandOddsRatio檢測(cè)是否有鏈偏差。
#MQRankSum：Mann-Whitney U test檢測(cè)是否有映射質(zhì)量的秩和差異。
#ReadPosRankSum：Mann-Whitney U test檢測(cè)是否有讀段位置的秩和差異。
#DP：每個(gè)變異位點(diǎn)的深度。
#--filter-name 'SNP_filter'：指定過(guò)濾標(biāo)簽的名稱為SNP_filter，這個(gè)標(biāo)簽會(huì)添加到vcf文件中的FILTER字段中，表示該記錄未通過(guò)過(guò)濾標(biāo)準(zhǔn)。
#--missing-values-evaluate-as-failing：指定如果某個(gè)變異位點(diǎn)缺少某個(gè)質(zhì)量指標(biāo)，則認(rèn)為該記錄未通過(guò)過(guò)濾標(biāo)準(zhǔn)。
#-O $WK/04Variants/02A1A2_SSHF_SNP_raw_filter.vcf：指定輸出的vcf文件。
這一部分使用SelectVariants功能，從02A1A2_SSHF_SNP_raw_filter.vcf文件中提取出通過(guò)過(guò)濾標(biāo)準(zhǔn)（即FILTER字段為PASS）的變異，并輸出到03A1A2_SSHF_pass_SNP.vcf文件中。

#03 提取過(guò)濾好的SNP 
/usr/bin/java -Xmx32G -jar $GATK/gatk-package-4.1.0.0-local.jar SelectVariants \ 
-R $REF/$fasta \
 -V $WK/04Variants/02A1A2_SSHF_SNP_raw_filter.vcf \ 
--exclude-filtered \ 
-O $WK/04Variants/03A1A2_SSHF_pass_SNP.vcf

每行代碼的含義如下：
#/usr/bin/java -Xmx32G：同上
#-jar $GATK/gatk-package-4.1.0.0-local.jar：同上。
#SelectVariants：同上。
#-R $REF/$fasta：同上。
#-V $WK/04Variants/02A1A2_SSHF_SNP_raw_filter.vcf：指定輸入的vcf文件，與上一步輸出相同。
#--exclude-filtered：指定排除掉未通過(guò)過(guò)濾標(biāo)準(zhǔn)（即FILTER字段不為PASS）的記錄。
#-O $WK/04Variants/03A1A2_SSHF_pass_SNP.vcf：指定輸出的vcf文件。
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容