首先我们先下载一组全外显子测序数据。nabi sra库,随机找了一个。
来自受试者“16177_CCPM_1300019”(SRR28391647, SRR28398576)的样本“16177_CCPM_1300019_BB5”的基因组DNA配对端文库“0369547849_Illumina_P5-Popal_P7-Hefel”的Illumina随机外显子测序
下载下来,转为两个配对的fq文件。过程可参考3.windows下Ubuntu,sratoolkit软件,从ncbi的sra数据库下载数据。_sratools下载数据-CSDN博客
这样我们得到了两个配对的fq文件,如果太大,可以压缩一下。
1.质控
我们直接进行fastp进行质量控制,这时我们可以实现以下流程化。
ls *_1.fq.gz >1
ls *_2.fq.gz >2
paste 1 2 > merge
sed -i "s/.fq.gz//g" merge
cat merge | while read id ; do
fastp -i "${id}_1.fq.gz" -I "${id}_2.fq.gz" -o fastp/"${id}_clean_1.fq.gz" -O fastp/"${id}_clean_2.fq.gz" -j fastp/"${id}.json" -h fastp/"${id}.html"
done
这样我们就可以对一个目录下的样本进行批量指控了,将它们填写到脚本中,可以直接运行。
2.比对并排序
接下来,我们使用bwa mem模式对经过质控的fq文件进行比对,生成bam文件,进行排序。
cat ../merge | while read id ; do
bwa mem -t 20 -R "@RG\tID:${id}\tLB:${id}\tPL:Illumina\tSM:${id}" /mnt/h/db/bwa.db/hg38.fa ./"${id}_clean_1.fq.gz" ./"${id}_clean_2.fq.gz" | samtools sort -@ 2 -o human/"${id}.bam"
done
这上面的-R参数很重要。-R 接的是 Read Group的字符串信息,这是一个非常重要的信息,以@RG开头,它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的,这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。在Read Group中,有如下几个信息非常重要:
(1) ID,这是Read Group的分组ID,一般设置为测序的lane ID(不同lane之间的测序过程认为是独立的),下机数据中我们都能看到这个信息的,一般都是包含在fastq的文件名中;
(2) PL,指的是所用的测序平台,这个信息不要随便写!特别是当我们需要使用GATK进行后续分析的时候,更是如此!这是一个很多新手都容易忽视的一个地方,在GATK中,PL只允许被设置为:ILLUMINA,SLX,SOLEXA,SOLID,454,LS454,COMPLETE,PACBIO,IONTORRENT,CAPILLARY,HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台,当然,如果实在不知道,那么必须设置为UNKNOWN,名字方面不区分大小写。如果你在分析的时候这里没设置正确,那么在后续使用GATK过程中可能会碰到类似如下的错误:
ERROR MESSAGE: The platform (xx) associated with read group GATKSAMReadGroupRecord @RG:xx is not a recognized platform.
这个时候你需要对比对文件的header信息进行重写,就会稍微比较麻烦。
我们上面的例子用的是PL:illumina。如果你的数据是CG测序的那么记得不要写成CG!而要写COMPLETE。
(3) SM,样本ID,同样非常重要,有时候我们测序的数据比较多的时候,那么可能会分成多个不同的lane分布测出来,这个时候SM名字就是可以用于区分这些样本。
(4) LB,测序文库的名字,这个重要性稍微低一些,主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到,如果上面的lane ID足够用于区分的话,也可以不用设置LB;
除了以上这四个之外,还可以自定义添加其他的信息,不过如无特殊的需要,对于序列比对而言,这4个就足够了。这些信息设置好之后,在RG字符串中要用制表符(\t)将它们分开。
3.标记PCR重复,使用picard
cat ../merge | while read id ; do
java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates \
I="${id}.bam" \
O="${id}.markup.bam" \
M="${id}.markdup_metrics.txt"
done
cat ../merge
: 这个命令会将merge
文件的内容输出到标准输出流。while read id; do ... done
: 这是一个while
循环,它会逐行读取cat
命令的输出,并将每一行的内容赋值给变量id
。java -jar /mnt/h/softwore/picard/picard.jar MarkDuplicates ...
: 这是 Picard 工具的命令,用于标记重复的 reads。具体的参数解释如下:-jar /mnt/h/softwore/picard/picard.jar
: 指定 Picard 工具的 JAR 文件路径。MarkDuplicates
: 执行标记重复 reads 的操作。I="${id}.bam"
: 输入 BAM 文件的路径和文件名,${id}.bam
表示根据当前循环的id
变量构建的 BAM 文件名。O="${id}.markup.bam"
: 输出标记重复 reads 后的 BAM 文件的路径和文件名,${id}.markup.bam
表示根据当前循环的id
变量构建的标记后的 BAM 文件名。M="${id}.markdup_metrics.txt"
: 输出标记重复 reads 后的统计信息文件的路径和文件名,${id}.markdup_metrics.txt
表示根据当前循环的id
变量构建的统计信息文件名。
4.samtools建立索引
cat ../merge | while read id ; do
samtools index "${id}.markup.bam" -@ 10
done
接下来,我们就可以进入GATK4的流程了。
这是我最近跑的一个流程,因为我在很久之前就跑过,所以这次没有很多基础信息,我稍后会分享一个资料,来对这些流程进行一个原理的解释。这个流程就是跑一下主要过程,想要建立一个流程,得要考虑很多东西,选择软件,修改添加参数啥的等,还是得有具体的项目实施。