GATA检测变异超详细教程

应该是最详细的生信分析教程,主要由新必应(也就是搭载的chat GP T 4.0)协助完成,我只提供了分析代码和做了最终调整排版。

ChatGPT对代码进行了逐行解释和分节进行总结说明,同时对我的代码进行了纠错(确实是我没有发现的小错误,挺惊喜,以后让AI进行改脚本应该效率挺高)

对于生信初学者绝对是很好的助攻,同样对生信分析流程记录和教程撰写来说也是全新的体验,省去了很多对分析流程的调整和解释说明,使用起来非常方便,特别是对不熟悉的代码的解读以及对代码进行纠错

虽然只是AI带来的编程学习和分析的一些便利,但总感觉这只是开始,后面不知道会给我们的学习和生活带来多少影响和颠覆


分析流程正文

1、定义变量和路径

第一部分定义了一些变量,表示你后面要用到的路径和文件


fasta=Malus_domestica_HFTH1_v1.0.a1.fa #定义参考基因组的文件名

TEMP=/data/temp #定义一个临时目录

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 #定义项目目录

2、Trimmomatic清洗数据

  • 第二部分检查是否有一个叫01data/clean的目录,如果没有就创建一个。这个目录会存放去除接头和低质量碱基后的修剪和清洗过的序列。
  • 使用Trimmomatic工具来对sampleid.txt中列出的每个样本的双端序列进行修剪和清洗。Trimmomatic需要一些参数来指定如何修剪序列,比如LEADING, TRAILING, SLIDINGWINDOW, MINLEN 和 TOPHRED33。
if [ ! -d ./01data/clean ] #如果没有./01data/clean这个目录

then mkdir -p ./01data/clean #就创建一个./01data/clean这个目录

fi

for i in cat $WK/01data/sampleid.txt #对$WK/01data/sampleid.txt文件中列出的每个样本名字做循环

#01rawdata

clean do trimmomatic PE -threads 12 \ #用trimmomatic工具对双端序列进行修剪和清洗,使用12个线程,并把结果输出到下面四个文件中,分别是配对和不配对的序列。

$WK/01data/${i}_1.fq.gz $WK/01data/${i}_2.fq.gz \ #输入原始双端序列文件,${i}表示样本名字。

$WK/01data/clean/${i}_1_paired_clean.fq.gz \ #输出修剪和清洗后的第一条配对序列文件。

$WK/01data/clean/${i}_1_unpair_clean.fq.gz \ #输出修剪和清洗后的第一条不配对序列文件。

$WK/01data/clean/${i}_2_paired_clean.fq.gz \ #输出修剪和清洗后的第二条配对序列文件。

$WK/01data/clean/${i}_2_unpair_clean.fq.gz \ #输出修剪和清洗后的第二条不配对序列文件。

LEADING:3 TRAILING:3 \ #指定去除每条序列开头和结尾质量值低于3的碱基。

SLIDINGWINDOW:4:20 MINLEN:50 TOPHRED33 #指定使用一个长度为4碱基、质量值阈值为20的滑动窗口来扫描每条序列,并去除窗口内平均质量值低于阈值的部分;指定去除长度小于50碱基的序列;指定使用Phred+33编码格式来表示质量值。

done

3、clean data数据质控

  • 检查是否有一个叫01data/clean/fastqc的目录,如果没有就创建一个。这个目录会存放修剪和清洗过的序列的质量控制报告。

  • 使用一个叫FastQC 的工具来为每对修剪和清洗过的序列生成质量控制报告。FastQC提供了一系列模块化的分析方法来检查序列质量的各个方面,比如每个碱基位置上的质量值、每条序列上的质量值、每个碱基位置上碱基含量、每个碱基位置上GC含量、每条序列上GC含量、每个碱基位置上N含量、序列长度分布、序列重复水平、过度表达序列等等。FastQC还生成HTML文件,包含图表和表格来可视化质量指标。

  • 使用MultiQC工具来把FastQC报告汇总成一个总结报告,可以方便地浏览。

if [ ! -d ./01data/clean/fastqc ] #如果没有./01clean/fastqc这个目录

then mkdir -p ./01clean/fastqc #就创建一个./clean/fastqc这个目录

fi

for i in cat $WK/01sampleid.txt #对$WK/sampleid.txt文件中列出的每个样本名字做循环

do

fastqc $WK/data/clean/${i}_1_paired_clean.fq.gz -o $WK/data/clean/fastqc/ #用fastqc工具对修剪和清洗后的第一条配对序列生成质量控制报告,并把结果输出到$WK/data/clean/fastqc/

fastqc $WK/data/clean/${i}_2_paired_clean.fq.gz -o $WK/data/clean/fastqc/ #用fastqc工具对修剪和清洗后的第二条配对序列生成质量控制报告,并把结果输出到$WK/data/clean/fastqc/

done

multiqc $WK/data/clean/fastqc/ -o $WK/data/clean/fastqc/ #用multiqc工具对fastqc生成的报告进行汇总,并把结果输出到$WK/data/clean/fastqc/

4、比对

  • 进行基因组比对的步骤,使用了bwa软件包来实现。

  • 并进行PCR重复序列的标记和删除的步骤,使用picard软件包来实现。

  • 然后进行比对结果的统计和评估的步骤,使用了samtools软件包来实现。

if [ ! -d $WK/02mapping ] # 如果不存在目录$WK/02mapping,就创建一个

then mkdir -p $WK/02mapping # 创建目录$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 # 对于每一个样本id,都执行以下操作

do

bwa mem -t 12 -M -R “@RG\tID:$id\tSM:$id\tLB:WES\tPL:Illumina” \ # 使用bwa mem命令将样本的双端测序数据与参考基因组进行比对,使用12个线程,标记次优比对位置,添加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文件,并将比对结果通过管道传递给samtools sort命令,使用12个线程将结果按照坐标排序并输出为bam格式文件,保存在$WK/02mapping/bwa_mem${id}.bam中java -Xmx64g -jar $PICARD MarkDuplicates \ # 使用java运行picard的MarkDuplicates工具,分配64g内存空间

I=$WK/02mapping/bwa_mem_${id}.bam \ # 指定输入的bam文件为上一步的比对结果文件

O=$WK/02mapping/${id}_sort_rmdup.bam \ # 指定输出的bam文件为去除重复后的文件,保存在$WK/02mapping/${id}_sort_rmdup.bam中CREATE_INDEX=true \ # 创建索引文件

REMOVE_DUPLICATES=true \ # 删除重复序列而不只是标记它们

M=$WK/02mapping/${id}_sort_rmdup_metrics.txt # 输出重复序列的统计信息到文本文件中,保存在$WK/02mapping/${id}_sort_rmdup_metrics.txt中samtools flagstat $WK/02mapping/${id}_sort_rmdup.bam >
$WK/02mapping/${id}_deduped_stat.txt # 使用samtools flagstat命令对去除重复后的bam文件进行统计分析,并将结果输出到文本文件中,保存在$WK/02mapping/${id}_deduped_stat.txt中samtools coverage $WK/02mapping/${id}_sort_rmdup.bam > $WK/02mapping/${id}_deduped_coverage.txt # 使用samtools coverage命令计算去除重复后的bam文件的覆盖度,并将结果输出到文本文件中,保存在$WK/02mapping/${id}_deduped_coverage.txt中

done

5、检测变异

  • 对每个样本进行HaplotypeCaller分析,并输出GVCF格式的文件。这个步骤可以检测出每个样本中的变异位点和非变异位点,并将非变异位点分成块以压缩文件大小。

  • 用CombineGVCFs模块将多个GVCF文件合并成一个大的GVCF文件,以便于后续进行联合分析。这个步骤可以将不同样本中的变异信息整合到一起,并保留每个样本的基因型似然值。

  • 用GenotypeGVCFs模块对合并后的GVCF文件进行联合基因型分析,输出最终的VCF文件。这个步骤可以根据所有样本的数据计算出每个位点上每个样本的最可能基因型,并给出质量评估和过滤标记。

#01检测变异

if [ ! -d $WK/03vcf ] # 如果不存在目录$WK/03vcf,就创建这个目录

then mkdir -p $WK/03vcf # 创建目录$WK/03vcf

fi # 结束if语句

for i in `cat $WK/01data/sampleid.txt` # 对$WK/01data/sampleid.txt文件中的每一行(即每个样本ID)执行循环

do # 开始循环

gatk --java-options "-Xmx32g -Djava.io.tmpdir=$TEMP" HaplotypeCaller \ # 用gatk软件包调用HaplotypeCaller模块,设置java选项为-Xmx32g(最大内存为为32G)和-Djava.io.tmpdir=$TEMP(临时文件夹为$TEMP)
     -R $REF/GDDH13_1-1_formatted2.fasta \ # 设置参考基因组文件为$REF/GDDH13_1-1_formatted2.fasta
     -I $WK/02mapping/${i}_sort_rmdup.bam \ # 设置输入文件为$WK/02mapping/${i}_sort_rmdup.bam,其中${i}是样本ID
     -O $WK/03vcf/${i}_gvcf.gz \ # 设置输出文件为$WK/03vcf/${i}_gvcf.gz,其中${i}是样本ID,输出格式为GVCF
     -ERC GVCF # 设置发射模式为GVCF,即只输出变异位点和非变异位点的块信息
 done # 结束循环

# 02合并多个GVCF文件
 ls $WK/03vcf/*_gvcf.gz > $WK/03vcf/gvcf_all.list # 列出所有以_gvcf.gz结尾的文件,并保存到$WK/03vcf/gvcf_all.list中
 /usr/bin/java -Xmx32G -jar $GATK/gatk-package-4.1.0.0-local.jar CombineGVCFs \ # 用/usr/bin/java软件调用$GATK/gatk-package-4.1.0.0-local.jar包中的CombineGVCFs模块,设置最大内存为32G 
     -R $REF/$fasta \ # 设置参考基因组文件为$REF/$fasta 
    
 -V $WK/02HaplotypeCaller_gvcf/bwa_mem_A1_1_g_vcf.gz \ # 设置输入文件之一为$WK/02HaplotypeCaller_gvcf/bwa_mem_A1_1_g_vcf.gz 
     -V $WK/02HaplotypeCaller_gvcf/bwa_mem_A2_gr_vcf.gz \ # 设置输入文件之一为$WK/02HaplotypeCaller_gvcf/bwa_mem_A2_gr_vcf.gz 
     -V $PK/parents_genome/ss_genome/03vcf/ss_gvcf.gz \  # 设置输入文件之一为$PK/parents_genome/ss_genome/03vcf/ss_gvcf.gz 
     -V $PK/parents_genome/hanfu1_genome/hanfu1_gvcf.gz \  # 设置输入文件之一为$PK /parents_genome/hanfu1_genome/hanfu1_gvcf.gz 
     -O $WK /03Combine_vcf/A1A2_SSHF_Combine_G.vcf  # 设置输出文件为$ 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软件调用$GATK/gatk-package-4.1.0.0-local.jar包中的GenotypeGVCFs模块,设置最大内存为32G和临时文件夹为$TEMP 
        -R $REF/$fasta \  # 设置参考基因组
        -V $WK/03Combine_vcf/A1A2_SSHF_Combine_G.vcf \ # 设置输入文件为$WK/03Combine_vcf/A1A2_SSHF_Combine_G.vcf,即合并后的GVCF文件
        -O $WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf  # 设置输出文件为$WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf,即最终的VCF文件

6、 提取SNP(InDel)并进行过滤

这一部分使用SelectVariants功能,从A1A2_SSHF_Genotype.vcf文件中提取出SNP类型的变异,并输出到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程序运行GATK,并设置最大内存为32GB。
# `-jar $GATK/gatk-package-4.1.0.0-local.jar`:指定运行的GATK版本为4.1.0.0。
# `SelectVariants`:指定使用SelectVariants功能。
# `-R $REF/$fasta`:指定参考基因组序列文件,$REF和$fasta是变量,需要根据实际情况替换。
# `-V $WK/03Combine_vcf/A1A2_SSHF_Genotype.vcf`:指定输入的vcf文件,$WK是变量,需要根据实际情况替换。
# `--select-type SNP`:指定提取SNP类型的变异。
# `-O $WK/04Variants/01A1A2_SSHF_SNP_raw.vcf`:指定输出的vcf文件。
这一部分使用VariantFiltration功能,对01A1A2_SSHF_SNP_raw.vcf文件中的SNP进行质量过滤,并添加过滤标签到02A1A2_SSHF_SNP_raw_filter.vcf文件中。

#02 过滤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`:指定过滤掉窗口大小为10bp内有3个或以上变异位点的记录。
# `--filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || DP < 8.0 :指定过滤表达式,根据一些质量指标对变异进行过滤。这些指标包括:
#QD:每个变异位点的质量值除以深度。
#MQ:比对到参考基因组的平均质量值。
#FS:Fisher’s exact test检测是否有分层偏差。
#SOR:StrandOddsRatio检测是否有链偏差。
#MQRankSum:Mann-Whitney U test检测是否有映射质量的秩和差异。
#ReadPosRankSum:Mann-Whitney U test检测是否有读段位置的秩和差异。
#DP:每个变异位点的深度。
#--filter-name 'SNP_filter':指定过滤标签的名称为SNP_filter,这个标签会添加到vcf文件中的FILTER字段中,表示该记录未通过过滤标准。
#--missing-values-evaluate-as-failing:指定如果某个变异位点缺少某个质量指标,则认为该记录未通过过滤标准。
#-O $WK/04Variants/02A1A2_SSHF_SNP_raw_filter.vcf:指定输出的vcf文件。
这一部分使用SelectVariants功能,从02A1A2_SSHF_SNP_raw_filter.vcf文件中提取出通过过滤标准(即FILTER字段为PASS)的变异,并输出到03A1A2_SSHF_pass_SNP.vcf文件中。

#03 提取过滤好的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:指定排除掉未通过过滤标准(即FILTER字段不为PASS)的记录。
#-O $WK/04Variants/03A1A2_SSHF_pass_SNP.vcf:指定输出的vcf文件。
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 216,324评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,356评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 162,328评论 0 353
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 58,147评论 1 292
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 67,160评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 51,115评论 1 296
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 40,025评论 3 417
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,867评论 0 274
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,307评论 1 310
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,528评论 2 332
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,688评论 1 348
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,409评论 5 343
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 41,001评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,657评论 0 22
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,811评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,685评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,573评论 2 353

推荐阅读更多精彩内容