生信軟件:fastp

摘錄自:

http://www.reibang.com/p/6f492058da5b
https://zhuanlan.zhihu.com/p/33601691

github:https://github.com/OpenGene/fastp

類似軟件:Trimmomatic(成熟淡诗,復(fù)雜)照藻,SOAPnuke(華大基因)

fastp的特性:

  • 對(duì)數(shù)據(jù)自動(dòng)進(jìn)行全方位質(zhì)控韵洋,生成人性化的報(bào)告
  • 過(guò)濾功能(低質(zhì)量,太短诱担,太多N……);
  • 對(duì)每一個(gè)序列的頭部或尾部,計(jì)算滑動(dòng)窗內(nèi)的質(zhì)量均值,并將均值較低的子序列進(jìn)行切除(類似Trimmomatic的做法,但是快非常多);
  • 全局剪裁 (在頭/尾部本股,不影響去重),對(duì)于Illumina下機(jī)數(shù)據(jù)往往最后一到兩個(gè)cycle需要這樣處理;
  • 去除接頭污染桐腌。厲害的是,你不用輸入接頭序列苟径,因?yàn)樗惴〞?huì)自動(dòng)識(shí)別接頭序列并進(jìn)行剪裁;
  • 對(duì)于雙端測(cè)序(PE)的數(shù)據(jù)案站,軟件會(huì)自動(dòng)查找每一對(duì)read的重疊區(qū)域,并對(duì)該重疊區(qū)域中不匹配的堿基對(duì)進(jìn)行校正;
  • 去除尾部的polyG棘街。對(duì)于Illumina NextSeq/NovaSeq的測(cè)序數(shù)據(jù)蟆盐,因?yàn)槭莾缮òl(fā)光,polyG是常有的事遭殉,所以該特性對(duì)該兩類測(cè)序平臺(tái)默認(rèn)打開;
  • 對(duì)于PE數(shù)據(jù)中的overlap區(qū)間中不一致的堿基對(duì)石挂,依據(jù)質(zhì)量值進(jìn)行校正;
  • 可以對(duì)帶分子標(biāo)簽(UMI)的數(shù)據(jù)進(jìn)行預(yù)處理,不管UMI在插入片段還是在index上险污,都可以輕松處理;
  • 可以將輸出進(jìn)行分拆痹愚,而且支持兩種模式,分別是指定分拆的個(gè)數(shù)蛔糯,或者分拆后每個(gè)文件的行數(shù);

以上功能大多都不需要輸入太多的參數(shù)拯腮,一些功能默認(rèn)已經(jīng)開啟,但是可以用參數(shù)關(guān)閉蚁飒。fastp完美支持gzip的輸入和輸出动壤,同時(shí)支持SE和PE數(shù)據(jù),而且不但支持像Illumina平臺(tái)的short read數(shù)據(jù)淮逻,也在一定程度上支持了PacBio/Nanopore的long reads數(shù)據(jù)琼懊。

fastp軟件會(huì)生成HTML格式的報(bào)告阁簸,而且該報(bào)告中沒有任何一張靜態(tài)圖片,所有的圖表都是使用JavaScript動(dòng)態(tài)繪制哼丈,非常具有交互性启妹。想要看一下樣板報(bào)告的,可以去以下鏈接:http://opengene.org/fastp/fastp.html

而且軟件的開發(fā)者還充分考慮到了各種自動(dòng)化分析的需求削祈,不但生成了人可讀的HTML報(bào)告翅溺,還生成了程序可讀性非常強(qiáng)的JSON結(jié)果,該JSON報(bào)告中的數(shù)據(jù)包含了HTML報(bào)告100%的信息髓抑,而且該JSON文件的格式還是特殊定制的咙崎,不但程序讀得爽,你用任何一款文本編輯器打開吨拍,一眼過(guò)去也會(huì)看得明明白白褪猛。想要看一下JSON結(jié)果長(zhǎng)什么樣的,可以去以下鏈接:http://opengene.org/fastp/fastp.json

下面我們先來(lái)看看fastp的具體參數(shù):

usage: fastp -i <in1> -o <out1> [-I <in1> -O <out2>] [options...]
options:
  # I/O options   即輸入輸出文件設(shè)置
  -i, --in1                          read1 input file name (string)
  -o, --out1                         read1 output file name (string [=])
  -I, --in2                          read2 input file name (string [=])
  -O, --out2                         read2 output file name (string [=])
  -6, --phred64                      indicates the input is using phred64 scoring (it'll be converted to phred33, so the output will still be phred33)
  -z, --compression                  compression level for gzip output (1 ~ 9). 1 is fastest, 9 is smallest, default is 2. (int [=2])
    --reads_to_process               specify how many reads/pairs to be processed. Default 0 means process all reads. (int [=0])
  
  # adapter trimming options   過(guò)濾序列接頭參數(shù)設(shè)置
  -A, --disable_adapter_trimming     adapter trimming is enabled by default. If this option is specified, adapter trimming is disabled
  -a, --adapter_sequence               the adapter for read1. For SE data, if not specified, the adapter will be auto-detected. For PE data, this is used if R1/R2 are found not overlapped. (string [=auto])
      --adapter_sequence_r2            the adapter for read2 (PE data only). This is used if R1/R2 are found not overlapped. If not specified, it will be the same as <adapter_sequence> (string [=])
    
  # global trimming options   剪除序列起始和末端的低質(zhì)量堿基數(shù)量參數(shù)
  -f, --trim_front1                  trimming how many bases in front for read1, default is 0 (int [=0])
  -t, --trim_tail1                   trimming how many bases in tail for read1, default is 0 (int [=0])
  -F, --trim_front2                  trimming how many bases in front for read2. If it's not specified, it will follow read1's settings (int [=0])
  -T, --trim_tail2                   trimming how many bases in tail for read2. If it's not specified, it will follow read1's settings (int [=0])

  # polyG tail trimming, useful for NextSeq/NovaSeq data   polyG剪裁
  -g, --trim_poly_g                  force polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data
      --poly_g_min_len                 the minimum length to detect polyG in the read tail. 10 by default. (int [=10])
  -G, --disable_trim_poly_g          disable polyG tail trimming, by default trimming is automatically enabled for Illumina NextSeq/NovaSeq data

  # polyX tail trimming
  -x, --trim_poly_x                    enable polyX trimming in 3' ends.
      --poly_x_min_len                 the minimum length to detect polyX in the read tail. 10 by default. (int [=10])
  
  # per read cutting by quality options   劃窗裁剪
  -5, --cut_by_quality5              enable per read cutting by quality in front (5'), default is disabled (WARNING: this will interfere deduplication for both PE/SE data)
  -3, --cut_by_quality3              enable per read cutting by quality in tail (3'), default is disabled (WARNING: this will interfere deduplication for SE data)
  -W, --cut_window_size              the size of the sliding window for sliding window trimming, default is 4 (int [=4])
  -M, --cut_mean_quality             the bases in the sliding window with mean quality below cutting_quality will be cut, default is Q20 (int [=20])
  
  # quality filtering options   根據(jù)堿基質(zhì)量來(lái)過(guò)濾序列
  -Q, --disable_quality_filtering    quality filtering is enabled by default. If this option is specified, quality filtering is disabled
  -q, --qualified_quality_phred      the quality value that a base is qualified. Default 15 means phred quality >=Q15 is qualified. (int [=15])
  -u, --unqualified_percent_limit    how many percents of bases are allowed to be unqualified (0~100). Default 40 means 40% (int [=40])
  -n, --n_base_limit                 if one read's number of N base is >n_base_limit, then this read/pair is discarded. Default is 5 (int [=5])
  
  # length filtering options   根據(jù)序列長(zhǎng)度來(lái)過(guò)濾序列
  -L, --disable_length_filtering     length filtering is enabled by default. If this option is specified, length filtering is disabled
  -l, --length_required              reads shorter than length_required will be discarded, default is 15. (int [=15])

  # low complexity filtering
  -y, --low_complexity_filter          enable low complexity filter. The complexity is defined as the percentage of base that is different from its next base (base[i] != base[i+1]).
  -Y, --complexity_threshold           the threshold for low complexity filter (0~100). Default is 30, which means 30% complexity is required. (int [=30])

  # filter reads with unwanted indexes (to remove possible contamination)
      --filter_by_index1               specify a file contains a list of barcodes of index1 to be filtered out, one barcode per line (string [=])
      --filter_by_index2               specify a file contains a list of barcodes of index2 to be filtered out, one barcode per line (string [=])
      --filter_by_index_threshold      the allowed difference of index barcode for index filtering, default 0 means completely identical. (int [=0])

  # base correction by overlap analysis options   通過(guò)overlap來(lái)校正堿基
  -c, --correction                   enable base correction in overlapped regions (only for PE data), default is disabled
  
  # UMI processing
  -U, --umi                          enable unique molecular identifer (UMI) preprocessing
      --umi_loc                      specify the location of UMI, can be (index1/index2/read1/read2/per_index/per_read, default is none (string [=])
      --umi_len                      if the UMI is in read1/read2, its length should be provided (int [=0])
      --umi_prefix                   if specified, an underline will be used to connect prefix and UMI (i.e. prefix=UMI, UMI=AATTCG, final=UMI_AATTCG). No prefix by default (string [=])
      --umi_skip                       if the UMI is in read1/read2, fastp can skip several bases following UMI, default is 0 (int [=0])

  # overrepresented sequence analysis
  -p, --overrepresentation_analysis    enable overrepresented sequence analysis.
  -P, --overrepresentation_sampling    One in (--overrepresentation_sampling) reads will be computed for overrepresentation analysis (1~10000), smaller is slower, default is 20. (int [=20])

  # reporting options
  -j, --json                         the json format report file name (string [=fastp.json])
  -h, --html                         the html format report file name (string [=fastp.html])
  -R, --report_title                 should be quoted with ' or ", default is "fastp report" (string [=fastp report])
  
  # threading options   設(shè)置線程數(shù)
  -w, --thread                       worker thread number, default is 3 (int [=3])
  
  # output splitting options
  -s, --split                        split output by limiting total split file number with this option (2~999), a sequential number prefix will be added to output name ( 0001.out.fq, 0002.out.fq...), disabled by default (int [=0])
  -S, --split_by_lines               split output by limiting lines of each file with this option(>=1000), a sequential number prefix will be added to output name ( 0001.out.fq, 0002.out.fq...), disabled by default (long [=0])
  -d, --split_prefix_digits          the digits for the sequential number padding (1~10), default is 4, so the filename will be padded as 0001.xxx, 0 to disable padding (int [=4])
  
  # help
  -?, --help                         print this message

雖然參數(shù)看起來(lái)比較多羹饰,但常用的主要包括以下幾個(gè)部分:

輸入輸出文件設(shè)置
接頭處理
全局裁剪(即直接剪掉起始和末端低質(zhì)量堿基)
滑窗質(zhì)量剪裁 (與trimmomatic相似)
過(guò)濾過(guò)短序列
校正堿基(用于雙端測(cè)序)
質(zhì)量過(guò)濾

功能

  1. 接頭處理
    fastp默認(rèn)啟用了接頭處理伊滋,但是可以使用-A命令來(lái)關(guān)掉。fastp可以自動(dòng)化地查找接頭序列并進(jìn)行剪裁队秩,也就是說(shuō)你可以不輸入任何的接頭序列笑旺,fastp全自動(dòng)搞定了!對(duì)于SE數(shù)據(jù)馍资,你還是可以-a參數(shù)來(lái)輸入你的接頭筒主,而對(duì)于PE數(shù)據(jù)則完全沒有必要,fastp基于PE數(shù)據(jù)的overlap分析可以更準(zhǔn)確地查找接頭鸟蟹,去得更干凈乌妙,而且對(duì)于一些接頭本身就有堿基不匹配情況處理得更好。fastp對(duì)于接頭去除會(huì)有一個(gè)匯總的報(bào)告建钥。

  2. 全局裁剪
    fastp可以對(duì)所有read在頭部和尾部進(jìn)行統(tǒng)一剪裁藤韵,該功能在去除一些測(cè)序質(zhì)量不好的cycle比較有用,比如151*2的PE測(cè)序中熊经,最后一個(gè)cycle通常質(zhì)量是非常低的泽艘,需要剪裁掉。使用-f和-t分別指定read1的頭部和尾部的剪裁镐依,使用-F和-T分別指定read2的頭部和尾部的剪裁悉盆。

  3. 滑窗質(zhì)量剪裁
    很多時(shí)候,一個(gè)read的低質(zhì)量序列都是集中在read的末端馋吗,也有少部分是在read的開頭焕盟。fastp支持像Trimmomatic那樣對(duì)滑動(dòng)窗口中的堿基計(jì)算平均質(zhì)量值,然后將不符合的滑窗直接剪裁掉。使用-5參數(shù)開啟在5’端脚翘,也就是read的開頭的剪裁灼卢,使用-3參數(shù)開啟在3’端,也就是read的末尾的剪裁来农。使用-W參數(shù)指定滑動(dòng)窗大小鞋真,默認(rèn)是4,使用-M參數(shù)指定要求的平均質(zhì)量值沃于,默認(rèn)是20涩咖,也就是Q20。

  4. 過(guò)濾過(guò)短序列
    默認(rèn)開啟多序列過(guò)濾繁莹,默認(rèn)值為15檩互,使用-L(--disable_length_filtering)禁止此默認(rèn)選項(xiàng)∽裳荩或使用-l(--length_required)自定義最短序列闸昨。

  5. 校正堿基(用于雙端測(cè)序)
    fastp支持對(duì)PE數(shù)據(jù)的每一對(duì)read進(jìn)行分析,查找它們的overlap區(qū)間薄风,然后對(duì)于overlap區(qū)間中不一致的堿基饵较,如果發(fā)現(xiàn)其中一個(gè)質(zhì)量非常高,而另一個(gè)非常低遭赂,則可以將非常低質(zhì)量的堿基改為相應(yīng)的非常高質(zhì)量值的堿基值循诉。此選項(xiàng)默認(rèn)關(guān)閉,可使用-c(--correction)開啟撇他。

  6. 質(zhì)量過(guò)濾
    fastp可以對(duì)低質(zhì)量序列茄猫,較多N的序列,該功能默認(rèn)是啟用的逆粹,但可以使用-Q參數(shù)關(guān)閉。使用-q參數(shù)來(lái)指定合格的phred質(zhì)量值炫惩,比如-q 15表示質(zhì)量值大于等于Q15的即為合格僻弹,然后使用-u參數(shù)來(lái)指定最多可以有多少百分比的質(zhì)量不合格堿基。比如-q 15 -u 40表示一個(gè)read最多只能有40%的堿基的質(zhì)量值低于Q15他嚷,否則會(huì)被扔掉蹋绽。使用-n可以限定一個(gè)read中最多能有多少個(gè)N。

例子
最后筋蓖,附一個(gè)簡(jiǎn)單的例子:

#!/bin/bash

for i in 74 75 76 82 83 84 85 86 87 88; do
    {
    fastp -i ~/RNAseq/cleandata/SRR17343${i}_1.fastq.gz -o SRR17343${i}_1.fastq.gz \
        -I ~/RNAseq/cleandata/SRR17343${i}_2.fastq.gz -O SRR17343${i}_2.fastq.gz \
        -Q --thread=5 --length_required=50 --n_base_limit=6 --compression=6
    }&
done
wait

fastp結(jié)果解讀

接下來(lái)卸耘,我們?cè)倏匆幌氯绾卫斫鈌astp生成的質(zhì)控報(bào)告。fastp的報(bào)告在單一文件中同時(shí)包含了過(guò)濾前和過(guò)濾后的統(tǒng)計(jì)結(jié)果粘咖,如果是PE數(shù)據(jù)蚣抗,則同時(shí)包含了read1和read2的統(tǒng)計(jì)結(jié)果。之前我們已經(jīng)說(shuō)過(guò)了瓮下,fastp會(huì)生成HTML的報(bào)告和JSON格式的報(bào)告翰铡。HTML報(bào)告的默認(rèn)文件名是fastp.html钝域,但是可以通過(guò)-h參數(shù)修改,JSON報(bào)告的默認(rèn)文件名是fastp.json锭魔,但是可以通過(guò)-j參數(shù)修改例证。而且fastp報(bào)告還有一個(gè)標(biāo)題,默認(rèn)是fastp report迷捧,這個(gè)也可以通過(guò)-R參數(shù)修改為你想要的標(biāo)題织咧。JSON格式的報(bào)告是優(yōu)化過(guò)的,人機(jī)皆可讀漠秋,適合進(jìn)階的用戶使用程序解析笙蒙,而這里我們重點(diǎn)關(guān)注HTML格式的報(bào)告。

  1. 質(zhì)量分布曲線圖
    我們第一關(guān)注的當(dāng)然是質(zhì)量膛堤,所以fastp提供了質(zhì)量分布曲線手趣,即每一個(gè)cycle的平均質(zhì)量值,而且fastp同時(shí)提供了A/T/C/G四種不同堿基的平均質(zhì)量肥荔,以及總的平均質(zhì)量绿渣,如下圖所示:


    image

    從上圖我們可以看到,一共有5條曲線燕耿,分別是A/T/C/G和mean中符。而且HTML報(bào)告中的每一個(gè)項(xiàng)目和分項(xiàng)目都是可以點(diǎn)擊進(jìn)行隱藏和展開的。

  2. 堿基含量分布曲線
    和質(zhì)量分布曲線類似誉帅,堿基含量分布曲線也是按照每一個(gè)cycle來(lái)的淀散,顯示了每一個(gè)位置的堿基含量。如下圖所示:


    image

    image

    從圖中可以看到蚜锨,fastp同時(shí)顯示了A/T/C/G/N/GC
    的每一個(gè)位置的比例和總的比例档插。而且如果你覺得頭部那里比較亂看不清的話,可以用鼠標(biāo)拉一個(gè)框亚再,它就放大了郭膛。

  3. KMER統(tǒng)計(jì)表格
    fastp對(duì)5個(gè)堿基長(zhǎng)度的所有組合的出現(xiàn)次數(shù)進(jìn)行了統(tǒng)計(jì),然后把它放在了一張表格中氛悬,表格的每一個(gè)元素為深背景白字则剃,背景越深,則表示重復(fù)次數(shù)越多如捅。這樣棍现,一眼望去,就可以發(fā)現(xiàn)有哪一些異常的信息镜遣。


    image

    從上面的KMER表格中己肮,我們可以發(fā)現(xiàn),GGGGG的顏色特別深,從鼠標(biāo)移上去之后顯示的信息中我們可以發(fā)現(xiàn)它的出現(xiàn)次數(shù)是平均次數(shù)的12.8倍朴肺,這是不正常的窖剑,因?yàn)镚GGGG的正常倍數(shù)應(yīng)該在1倍左右。幸好我們有fastp戈稿,它可以過(guò)濾掉這種polyG西土,讓數(shù)值較多地回歸正常。

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末鞍盗,一起剝皮案震驚了整個(gè)濱河市需了,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌般甲,老刑警劉巖肋乍,帶你破解...
    沈念sama閱讀 216,372評(píng)論 6 498
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場(chǎng)離奇詭異敷存,居然都是意外死亡墓造,警方通過(guò)查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,368評(píng)論 3 392
  • 文/潘曉璐 我一進(jìn)店門锚烦,熙熙樓的掌柜王于貴愁眉苦臉地迎上來(lái)觅闽,“玉大人,你說(shuō)我怎么就攤上這事涮俄◎茸荆” “怎么了?”我有些...
    開封第一講書人閱讀 162,415評(píng)論 0 353
  • 文/不壞的土叔 我叫張陵彻亲,是天一觀的道長(zhǎng)孕锄。 經(jīng)常有香客問(wèn)我,道長(zhǎng)苞尝,這世上最難降的妖魔是什么畸肆? 我笑而不...
    開封第一講書人閱讀 58,157評(píng)論 1 292
  • 正文 為了忘掉前任,我火速辦了婚禮宙址,結(jié)果婚禮上轴脐,老公的妹妹穿的比我還像新娘。我一直安慰自己曼氛,他們只是感情好豁辉,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,171評(píng)論 6 388
  • 文/花漫 我一把揭開白布令野。 她就那樣靜靜地躺著舀患,像睡著了一般。 火紅的嫁衣襯著肌膚如雪气破。 梳的紋絲不亂的頭發(fā)上聊浅,一...
    開封第一講書人閱讀 51,125評(píng)論 1 297
  • 那天,我揣著相機(jī)與錄音,去河邊找鬼低匙。 笑死旷痕,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的顽冶。 我是一名探鬼主播欺抗,決...
    沈念sama閱讀 40,028評(píng)論 3 417
  • 文/蒼蘭香墨 我猛地睜開眼,長(zhǎng)吁一口氣:“原來(lái)是場(chǎng)噩夢(mèng)啊……” “哼强重!你這毒婦竟也來(lái)了绞呈?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 38,887評(píng)論 0 274
  • 序言:老撾萬(wàn)榮一對(duì)情侶失蹤间景,失蹤者是張志新(化名)和其女友劉穎佃声,沒想到半個(gè)月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體倘要,經(jīng)...
    沈念sama閱讀 45,310評(píng)論 1 310
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡圾亏,尸身上長(zhǎng)有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,533評(píng)論 2 332
  • 正文 我和宋清朗相戀三年,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了封拧。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片志鹃。...
    茶點(diǎn)故事閱讀 39,690評(píng)論 1 348
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡,死狀恐怖哮缺,靈堂內(nèi)的尸體忽然破棺而出弄跌,到底是詐尸還是另有隱情,我是刑警寧澤尝苇,帶...
    沈念sama閱讀 35,411評(píng)論 5 343
  • 正文 年R本政府宣布铛只,位于F島的核電站,受9級(jí)特大地震影響糠溜,放射性物質(zhì)發(fā)生泄漏淳玩。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,004評(píng)論 3 325
  • 文/蒙蒙 一非竿、第九天 我趴在偏房一處隱蔽的房頂上張望蜕着。 院中可真熱鬧,春花似錦红柱、人聲如沸承匣。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,659評(píng)論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽(yáng)韧骗。三九已至,卻和暖如春零聚,著一層夾襖步出監(jiān)牢的瞬間袍暴,已是汗流浹背些侍。 一陣腳步聲響...
    開封第一講書人閱讀 32,812評(píng)論 1 268
  • 我被黑心中介騙來(lái)泰國(guó)打工, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留政模,地道東北人岗宣。 一個(gè)月前我還...
    沈念sama閱讀 47,693評(píng)論 2 368
  • 正文 我出身青樓,卻偏偏與公主長(zhǎng)得像淋样,于是被迫代替她去往敵國(guó)和親耗式。 傳聞我的和親對(duì)象是個(gè)殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,577評(píng)論 2 353

推薦閱讀更多精彩內(nèi)容