背景

前面介绍了如何从sra文件得到fastq文件

[Linux|生信]project4_01:批量下载sra文件并转化为fastq文件

今天介绍一下如何对fastq文件进行批量质控。在此申明，该系列处理的是单细胞转录组的数据，不是常规的Bulk转录组，故质控后的处理思路有所不同。

质控（QC）

# 一：质控前的初步看测序数据质量：fastqc与multiqc
# 1.激活前面生成的虚拟环境，进行fastqc与multiqc
conda activate biotree 
# 2.先进行fastqc, 每个*.fq.gz文件都会得到一个网页文件，可直接通过浏览器查看测序的质量
nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log &  
# 3.对fastqc后的zip数据进行multiqc，将多个质控报告合并为一个，便于查看
nohup multiqc ./*.zip -o ./ > ./multiqc.log &

fastqc结果图

注：

1. 需要注意“Encoding”的值，后续使用“trimmgalore”进行质控时有个参数根据该值确定
2. fastaqc报告解读，可以参考链接 multiqc结果图

去除低质量碱基

ls *gz | xargs -n 2 -t nohup trim_galore  -q 25 --phred33 --length 36 --stringency 3 --paired -j 4 -o ./  &

# 其他学员反馈如下代码也可行（未验证）
nohup trim_galore -q 25 --phred33 --stringency 3 --length 36 --paired *fastq.gz --o ./ &

命令及参数说明

上述命令需分为四部分理解：

1. ls *gz 星号为通配符，意思为匹配多个任意字符，该部分为列出当前目录下gz结尾的所有文件
2. | 竖线为管道符，作用为将前面的标准输出通过管道符作为后面命令的标准输入
3. xargs -n 2 -t xargs作用为定义前面的标准输出如何应用于后面(其他使用方法见文末链接) -n 为指定每次将2个文件名作为参数供后续使用，因为后续命令为处理双端序列，输入文件有2个 -t 作用为将待执行的命令显示于终端，可以判断命令是否为我们预想的命令
4. nohup trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -j 4 -o ./ & 其中nohup & 的作用前面的推文介绍过，是将中间的命令“trim_galore -q 25 --phred33 --length 36 --stringency 3 --paired -j 4 -o ./ ”置于后台执行，防止意外中断。 trim_galore 首先去除低质量碱基，然后去除3' 末端的adapter(如果没有指定具体的adapter，程序会自动检测前1million的序列） -q 设定Phred quality score阈值，默认为Phred 20 切除质量得分低于设定值的序列 --phred33 使用ASCII+33/64质量得分作为Phred得分选择-phred33或者-phred64，表示-测序平台使用的Phred quality score。需要确认:Sanger/Illumina 1.9+ encoding为33; Illumina 1.5 encoding为64。（从fastaqc质控报告“Basic Statistics”一栏中可知“Encoding”参数） --length 长度低于指定阈值的reads被丢掉 -s/--stringency 设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到 --paired 双端序列必须指定的参数 -j 指定并行的线程数，可加速任务 -o 后接输出文件

• 结果发现经过“trim_galore”去除低质量碱基后和原始fastq.gz的质控报告差别不大。主要在“Per Base Sequence Content”、“Sequence Duplication Levels”、“Per Sequence GC Content”上大部分都是“failed”。对于“Per Base Sequence Content”我准备利用“trimmomatic”去除前10个不好的碱基，对于“Sequence Duplication Levels”中的大量重复我查阅资料发现在单细胞数据中是正常的，对于最后GC含量异常那个指标还没什么处理思路。

• 随后参阅资料发现单细胞转录组与Bulk 转录组对质控报告提示的问题的处理方式是不同的。
• 单细胞转录组是对少量细胞进行的测序，是细胞层面表达量的衡量工具，更能反映细胞的异质性。Bulk 转录组是样本层面，对所有细胞取平均之后的表达量。
• “Sequence Duplication Levels”上表现出的重复很有可能是相关基因的高表达，是有意义的，不能贸然去除。在Biostars上也有广泛的讨论，总结为对单细胞数据不建议去除。
• 对于“Per Base Sequence Content”表现出来reads头部10个bp ,ATGC比例不平衡，如果是普通的转录组数据，可以采用下述“trimmomatic”就行去除，在此次的单细胞数据中不进行处理。

java -jar $HOME/Download/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 8 -phred33 SRR11618624_1_val_1.fq.gz SRR11618624_2_val_2.fq.gz -baseout test_clean.fq.gz HEADCROP:10

注：HEADCROP后以冒号接参数，不是常规的空格，有空格会报错

• 对于“Per Sequence GC Content”表示出的异常，暂时没有解决思路，大家有什么建议可以在评论区讨论。

• 最后对于普通（Bulk）转录组质控指标不理想的相应处理措施，网上的教程大都没有汇总。在此找到Fastqc官网给出的应对措施，供感兴趣的朋友研究。我是一个链接，点我。

• 推文多平台同步发布，公众号内容食用更佳

• 更多内容，请关注微信公众号“生信矿工”

• 如有意见或建议可以在评论区讨论

参考链接

1. xargs命令详细教程