转录组重建系统发育(一)使用fastp多rawdata进行处理

rawdata可以是从公司直接获得的测序数据或者是网上下载的数据,网上下载数据的流程参考使用SRA Toolkit下载NCBI-SRA原始数据
常规的质控和过滤数据是fastqc+trimmomatic,新出的fastp速度更快,操作更简单,而且一次完成质控过滤和出图。fastp是发表在《bioinfomatics》上的测序质控软件,文章链接:https://academic.oup.com/bioinformatics/article/34/17/i884/5093234,可以直接通过扫描fastq文件完成从质控到处理的工作,详细教程在:https://github.com/OpenGene/fastp

对数据自动进行全方位质控,生成人性化的报告。
过滤功能(低质量,太短,太多N……)。
对每一个序列的头部或尾部,计算滑动窗内的质量均值,并将均值较低的子序列进行切除(类似Trimmomatic的做法,但是快非常多)。
全局剪裁 (在头/尾部,不影响去重),对于Illumina下机数据往往最后一到两个cycle需要这样处理。
去除接头污染。厉害的是,你不用输入接头序列,因为算法会自动识别接头序列并进行剪裁。
对于双端测序(PE)的数据,软件会自动查找每一对read的重叠区域,并对该重叠区域中不匹配的碱基对进行校正。
去除尾部的polyG。对于Illumina NextSeq/NovaSeq的测序数据,因为是两色法发光,polyG是常有的事,所以该特性对该两类测序平台默认打开。
对于PE数据中的overlap区间中不一致的碱基对,依据质量值进行校正
可以对带分子标签(UMI)的数据进行预处理,不管UMI在插入片段还是在index上,都可以轻松处理。
可以将输出进行分拆,而且支持两种模式,分别是指定分拆的个数,或者分拆后每个文件的行数。

1. 软件安装

fastp安装和使用都非常简单,有conda的可以直接使用

conda install -c bioconda fastp

进行安装
也可以自己下载安装

wget http://opengene.org/fastp/fastp
chmod a+x ./fastp

2. 软件使用

fastp使用非常简单
单端测序输入命令

fastp -i SRS7102712.fastq.gz -o SRS7102712.clean.fastq.gz
#-i 后接输入文件名 -o 后接输出文件名

双端命令

fastp -i SRS7102712_1.fastq.gz -I SRS7102712_2.fastq.gz -o SRS7102712_1.clean.fastq.gz -O SRS7102712_2.clean.fastq.gz
#注意大小写 i和o的大小写分别代表两个输入文件和输出文件

3. 功能介绍和结果解读

获得如下图所示的4个文件

image.png

html后缀的文件可以使用图文界面查看,样板信息可以查看http://opengene.org/fastp/fast
获得的新的fastq.gz文件就是cleandata,使用这个就可以直接进行转录组组装了。

3.1过滤

fastp可以对低质量序列,较多N的序列,该功能默认是启用的,但可以使用-Q参数关闭。使用-q参数来指定合格的phred质量值,比如-q 15表示质量值大于等于Q15的即为合格,然后使用-u参数来指定最多可以有多少百分比的质量不合格碱基。比如-q 15 -u 40表示一个read最多只能有40%的碱基的质量值低于Q15,否则会被扔掉。使用-n可以限定一个read中最多能有多少个N。
fastp还默认启用了read长度过滤,但可以使用-L参数关闭。使用-l参数指定最低要求一个read有多长,比如-l 30表示低于30个碱基的read会被扔掉。这个功能可以用于实现常用的discard模式,以保证所有输出的序列都一样长

在fastp的HTML报告中,最头上的Summary表格很清楚地显示了过滤的统计信息,如下图所示:


image.png

3.2 接头处理

接头(adapter)污染的处理是FASTQ文件预处理中很重要的一步。fastp默认启用了接头处理,但是可以使用-A命令来关掉。fastp可以自动化地查找接头序列并进行剪裁,也就是说你可以不输入任何的接头序列,fastp全自动搞定了!对于SE数据,你还是可以-a参数来输入你的接头,而对于PE数据则完全没有必要,fastp基于PE数据的overlap分析可以更准确地查找接头,去得更干净,而且对于一些接头本身就有碱基不匹配情况处理得更好。fastp对于接头去除会有一个汇总的报告,如下图所示:

image.png

5.3滑窗质量剪裁

很多时候,一个read的低质量序列都是集中在read的末端,也有少部分是在read的开头。fastp支持像Trimmomatic那样对滑动窗口中的碱基计算平均质量值,然后将不符合的滑窗直接剪裁掉。使用-5参数开启在5’端,也就是read的开头的剪裁,使用-3参数开启在3’端,也就是read的末尾的剪裁。使用-W参数指定滑动窗大小,默认是4,使用-M参数指定要求的平均质量值,默认是20,也就是Q20。

5.4 PE数据的碱基校正

fastp支持对PE数据的每一对read进行分析,查找它们的overlap区间,然后对于overlap区间中不一致的碱基,如果发现其中一个质量非常高,而另一个非常低,则可以将非常低质量的碱基改为相应的非常高质量值的碱基值,如下图所示:


image.png

上图中所示的标红的T碱基是低质量序列,和高质量的A不匹配,它会被校正为A。该校正功能默认没有开启使用-c参数可以启用,对于一些对噪声容忍度低的应用,比如液体活检,建议开启。

5.5全局剪裁

fastp可以对所有read在头部和尾部进行统一剪裁,该功能在去除一些测序质量不好的cycle比较有用,比如151*2的PE测序中,最后一个cycle通常质量是非常低的,需要剪裁掉。使用-f和-t分别指定read1的头部和尾部的剪裁,使用-F和-T分别指定read2的头部和尾部的剪裁。

5.6 polyG剪裁

对于两色发光法的Illumina设备(NextSeq /
NovaSeq),因为在没有光信号情况下base calling的结果会返回G,所以在序列的尾端可能会出现较多的polyG,需要被去除。fastp会自动化地识别NextSeq / NovaSeq的数据,然后进行polyG识别和剪裁。如果你想强制开启该功能,可以指定-g参数,如果想强制关闭该功能,则可以指定-G参数。

5.7分子标签UMI处理

UMI在处理ctDNA类似的超低频突变检测应用中是十分有用的,为了更好地对带UMI的FASTQ文件进行预处理,fastp也很好地支持了UMI预处理功能。该功能默认没有启用,需要使用-U参数开启,另外需要使用--umi_loc来指定UMI所在的位置,它可以是(index1、 index2、 read1、 read2、 per_index、 per_read )中的一种,分别表示UMI是在index位置上,还是在插入片段中。如果指定了是在插入序列中,还需要使用 --umi_len 参数来指定UMI所占的碱基长度。

5.8输出文件切分

很多时候我们需要对输出的FASTQ进行切分,分成大小均匀的多个文件,这样可以使用比对软件并行地比对,提高并行处理的速度。fastp软件也提供了相应的功能,并且支持了两种模式,分别是使用参数-s指定切分后文件的个数,或者使用-S参数指定每个切分后文件的行数。

6质控报告解读

接下来,我们再看一下如何理解fastp生成的质控报告。fastp的报告在单一文件中同时包含了过滤前和过滤后的统计结果,如果是PE数据,则同时包含了read1和read2的统计结果。之前我们已经说过了,fastp会生成HTML的报告和JSON格式的报告。HTML报告的默认文件名是fastp.html,但是可以通过-h参数修改,JSON报告的默认文件名是fastp.json,但是可以通过-j参数修改。而且fastp报告还有一个标题,默认是fastp report,这个也可以通过-R参数修改为你想要的标题。JSON格式的报告是优化过的,人机皆可读,适合进阶的用户使用程序解析,而这里我们重点关注HTML格式的报告。

6.1质量分布曲线图

我们第一关注的当然是质量,所以fastp提供了质量分布曲线,即每一个cycle的平均质量值,而且fastp同时提供了A/T/C/G四种不同碱基的平均质量,以及总的平均质量,如下图所示:


image.png

从上图我们可以看到,一共有5条曲线,分别是A/T/C/G和mean。而且HTML报告中的每一个项目和分项目都是可以点击进行隐藏和展开的。

6.2碱基含量分布曲线

和质量分布曲线类似,碱基含量分布曲线也是按照每一个cycle来的,显示了每一个位置的碱基含量。如下图所示:


image

image

从图中可以看到,fastp同时显示了A/T/C/G/N/GC
的每一个位置的比例和总的比例。而且如果你觉得头部那里比较乱看不清的话,可以用鼠标拉一个框,它就放大了。

6.3 KMER统计表格

fastp对5个碱基长度的所有组合的出现次数进行了统计,然后把它放在了一张表格中,表格的每一个元素为深背景白字,背景越深,则表示重复次数越多。这样,一眼望去,就可以发现有哪一些异常的信息。
image

从上面的KMER表格中,我们可以发现,GGGGG的颜色特别深,从鼠标移上去之后显示的信息中我们可以发现它的出现次数是平均次数的12.8倍,这是不正常的,因为GGGGG的正常倍数应该在1倍左右。幸好我们有fastp,它可以过滤掉这种polyG,让数值较多地回归正常。

6.4过表达序列(overrepresented sequence)

fastp的最新版本(v0.12)提供了overrepresetned

sequence的分析,而且不但提供了这些overrepresented sequence的序列个数和占比,还提供了他们在测序cycles中的分布情况,这有利于分析各种问题。具体示例如下图所示:
image.png

fastp软件还在不断更新中,目前每星期都有新功能开发出来,所以要想了解fastp软件的最新动态,请关注该软件的github项目地址https://github.com/OpenGene/fas
知乎参考链接:https://zhuanlan.zhihu.com/p/33601691

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 200,045评论 5 468
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 84,114评论 2 377
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 147,120评论 0 332
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 53,902评论 1 272
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 62,828评论 5 360
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,132评论 1 277
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,590评论 3 390
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,258评论 0 254
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,408评论 1 294
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,335评论 2 317
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,385评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,068评论 3 315
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,660评论 3 303
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,747评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 30,967评论 1 255
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,406评论 2 346
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 41,970评论 2 341

推荐阅读更多精彩内容