应该是最详细的生信分析教程,主要由新必应(也就是搭载的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文件。