windows ubuntu 子系统：肿瘤全外篇，2. fq 数据质控，比对。

news2025/1/11 0:50:55

1.质控

2.比对并排序

3.标记PCR重复,使用picard

4.samtools建立索引

首先我们先下载一组全外显子测序数据。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 变量构建的统计信息文件名。