J

Jay2Echo

V1

2022/12/07阅读:57主题:默认主题

[Linux|生信]project4_02:质控过滤

背景

前面介绍了如何从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结果图

fastqc结果图 注:

  1. 需要注意“Encoding”的值,后续使用“trimmgalore”进行质控时有个参数根据该值确定
  2. fastaqc报告解读,可以参考链接 multiqc结果图
multiqc结果图
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含量异常那个指标还没什么处理思路。
“Per Base Sequence Content” Failed
“Per Base Sequence Content” Failed
“Per Sequence GC Content” mostly Failed
“Per Sequence GC Content” mostly Failed
“Sequence Duplication Levels” Failed
“Sequence Duplication Levels” Failed
  • 随后参阅资料发现单细胞转录组与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命令详细教程

分类:

后端

标签:

大数据

作者介绍

J
Jay2Echo
V1