轉錄組數(shù)據(jù)分析—fastp v0.23.1

用于執(zhí)行FASTQ數(shù)據(jù) 的質量控制唐全,支持單端和雙端短讀數(shù)據(jù),作為一體化的 FASTQ預處理器垃帅,fastp提供了 諸如質量分析剪勿,接頭修整酱固,讀 取過濾和基本校正等功能头朱。 在處理完所有讀取后班眯,這些值 將合并署隘,并且報告者將生成 HTML和JSON格式的報告 亚隙。它可以僅僅掃描FASTQ文件一次诊霹,就完成比FASTQC+cutadapt+Trimmomatic這三個軟件加起來的功能還多很多的功能畅哑,而且速度上比僅僅使用Trimmomatic一個軟件還要快3倍左右荠呐,因為它使用C++開發(fā)泥张,處處使用了高效算法媚创,而且完美支持多線程钞钙!


Fig.1
Fig.2

安裝

使用conda安裝conda install -c bioconda fastp=0.23.2瘫怜。

#看一下是否成功鲸湃,順便看一下help文檔
fastp
fastp: an ultra-fast all-in-one FASTQ preprocessor
version 0.23.2
usage: fastp [options] ...
options:
  -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 [=])
      --unpaired1                      for PE input, if read1 passed QC but read2 not, it will be written to unpaired1\. Default is to discard it. (string [=])
      --unpaired2                      for PE input, if read2 passed QC but read1 not, it will be written to unpaired2\. If --unpaired2 is same as --unpaired1 (default mode), both unpaired reads will be written to this same file. (string [=])
      --overlapped_out                 for each read pair, output the overlapped region if it has no any mismatched base. (string [=])
      --failed_out                     specify the file to store reads that cannot pass the filters. (string [=])
  -m, --merge                          for paired-end input, merge each pair of reads into a single read if they are overlapped. The merged reads will be written to the file given by --merged_out, the unmerged reads will be written to the files specified by --out1 and --out2\. The merging mode is disabled by default.
      --merged_out                     in the merging mode, specify the file name to store merged output, or specify --stdout to stream the merged output (string [=])
      --include_unmerged               in the merging mode, write the unmerged or unpaired reads to the file specified by --merge. Disabled by default.
  -6, --phred64                        indicate 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 4\. (int [=4])
      --stdin                          input from STDIN. If the STDIN is interleaved paired-end FASTQ, please also add --interleaved_in.
      --stdout                         stream passing-filters reads to STDOUT. This option will result in interleaved FASTQ output for paired-end output. Disabled by default.
      --interleaved_in                 indicate that <in1> is an interleaved FASTQ which contains both read1 and read2\. Disabled by default.
      --reads_to_process               specify how many reads/pairs to be processed. Default 0 means process all reads. (int [=0])
      --dont_overwrite                 don't overwrite existing files. Overwritting is allowed by default.
      --fix_mgi_id                     the MGI FASTQ ID format is not compatible with many BAM operation tools, enable this option to fix it.
  -V, --verbose                        output verbose log information (i.e. when every 1M reads are processed).
  -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 [=auto])
      --adapter_fasta                  specify a FASTA file to trim both read1 and read2 (if PE) by all the sequences in this FASTA file (string [=])
      --detect_adapter_for_pe          by default, the auto-detection for adapter is for SE data input only, turn on this option to enable it for PE data.
  -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])
  -b, --max_len1                       if read1 is longer than max_len1, then trim read1 at its tail to make it as long as max_len1\. Default 0 means no limitation (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])
  -B, --max_len2                       if read2 is longer than max_len2, then trim read2 at its tail to make it as long as max_len2\. Default 0 means no limitation. If it's not specified, it will follow read1's settings (int [=0])
  -D, --dedup                          enable deduplication to drop the duplicated reads/pairs
      --dup_calc_accuracy              accuracy level to calculate duplication (1~6), higher level uses more memory (1G, 2G, 4G, 8G, 16G, 24G). Default 1 for no-dedup mode, and 3 for dedup mode. (int [=0])
      --dont_eval_duplication          don't evaluate duplication rate to save time and use less memory.
  -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
  -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])
  -5, --cut_front                      move a sliding window from front (5') to tail, drop the bases in the window if its mean quality < threshold, stop otherwise.
  -3, --cut_tail                       move a sliding window from tail (3') to front, drop the bases in the window if its mean quality < threshold, stop otherwise.
  -r, --cut_right                      move a sliding window from front to tail, if meet one window with mean quality < threshold, drop the bases in the window and the right part, and then stop.
  -W, --cut_window_size                the window size option shared by cut_front, cut_tail or cut_sliding. Range: 1~1000, default: 4 (int [=4])
  -M, --cut_mean_quality               the mean quality requirement option shared by cut_front, cut_tail or cut_sliding. Range: 1~36 default: 20 (Q20) (int [=20])
      --cut_front_window_size          the window size option of cut_front, default to cut_window_size if not specified (int [=4])
      --cut_front_mean_quality         the mean quality requirement option for cut_front, default to cut_mean_quality if not specified (int [=20])
      --cut_tail_window_size           the window size option of cut_tail, default to cut_window_size if not specified (int [=4])
      --cut_tail_mean_quality          the mean quality requirement option for cut_tail, default to cut_mean_quality if not specified (int [=20])
      --cut_right_window_size          the window size option of cut_right, default to cut_window_size if not specified (int [=4])
      --cut_right_mean_quality         the mean quality requirement option for cut_right, default to cut_mean_quality if not specified (int [=20])
  -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])
  -e, --average_qual                   if one read's average quality score <avg_qual, then this read/pair is discarded. Default 0 means no requirement (int [=0])
  -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])
      --length_limit                   reads longer than length_limit will be discarded, default 0 means no limitation. (int [=0])
  -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_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])
  -c, --correction                     enable base correction in overlapped regions (only for PE data), default is disabled
      --overlap_len_require            the minimum length to detect overlapped region of PE reads. This will affect overlap analysis based PE merge, adapter trimming and correction. 30 by default. (int [=30])
      --overlap_diff_limit             the maximum number of mismatched bases to detect overlapped region of PE reads. This will affect overlap analysis based PE merge, adapter trimming and correction. 5 by default. (int [=5])
      --overlap_diff_percent_limit     the maximum percentage of mismatched bases to detect overlapped region of PE reads. This will affect overlap analysis based PE merge, adapter trimming and correction. Default 20 means 20%. (int [=20])
  -U, --umi                            enable unique molecular identifier (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])
  -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])
  -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])
  -w, --thread                         worker thread number, default is 3 (int [=3])
  -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])
      --cut_by_quality5                DEPRECATED, use --cut_front instead.
      --cut_by_quality3                DEPRECATED, use --cut_tail instead.
      --cut_by_quality_aggressive      DEPRECATED, use --cut_right instead.
      --discard_unmerged               DEPRECATED, no effect now, see the introduction for merging.
  -?, --help                           print this message

功能

  • 對數(shù)據(jù)自動進行全方位質控听系,生成人性化的報告靠胜。

  • 過濾功能(低質量毕源,太短霎褐,太多N……)冻璃。

  • 對每一個序列的頭部或尾部,計算滑動窗內的質量均值嫁审,并將均值較低的子序列進行切除(類似Trimmomatic的做法赖晶,但是快非常多)捂贿。

  • 全局剪裁 (在頭/尾部厂僧,不影響去重),對于Illumina下機數(shù)據(jù)往往最后一到兩個cycle需要這樣處理德召。

  • 去除接頭污染福荸。厲害的是敬锐,你不用輸入接頭序列台夺,因為算法會自動識別接頭序列并進行剪裁颤介。

  • 對于雙端測序(PE)的數(shù)據(jù)滚朵,軟件會自動查找每一對read的重疊區(qū)域辕近,并對該重疊區(qū)域中不匹配的堿基對進行校正移宅。

  • 去除尾部的polyG吞杭。對于Illumina NextSeq/NovaSeq的測序數(shù)據(jù),因為是兩色法發(fā)光绢掰,polyG是常有的事滴劲,所以該特性對該兩類測序平臺默認打開。

  • fastp支持對PE數(shù)據(jù)的每一對read進行分析萧芙,查找它們的overlap區(qū)間双揪,然后對于overlap區(qū)間中不一致的堿基渔期,如果發(fā)現(xiàn)其中一個質量非常高疯趟,而另一個非常低信峻,則可以將非常低質量的堿基改為相應的非常高質量值的堿基值,該校正功能默認沒有開啟使用-c參數(shù)可以啟用矾策,對于一些對噪聲容忍度低的應用贾虽,比如液體活檢蓬豁,建議開啟取募。

  • 可以對帶分子標簽(UMI)的數(shù)據(jù)進行預處理蟆技,不管UMI在插入片段還是在index上玩敏,都可以輕松處理。

  • 可以將輸出進行分拆质礼,而且支持兩種模式旺聚,分別是指定分拆的個數(shù),或者分拆后每個文件的行數(shù)眶蕉。

使用

fastp使用非常簡單砰粹,如果不清楚具體參數(shù)該如何設置,直接定義好輸入及輸出數(shù)據(jù)就可以快速實現(xiàn)造挽。

##單端數(shù)據(jù)
fastp -i R1.fastq.gz -o R1_clean.fastq.gz
##雙端數(shù)據(jù)
fastp -i R1.fastq.gz -o R1_clean.fastq.gz -I R2.fastq.gz -O R2_clean.fastq.gz

本次測試使用之前下載好的原始數(shù)據(jù)(GSE176393)碱璃,我的代碼

#批量操作
cat SRR_Acc_List.txt | while read line;
do
fastp -i 0-rawdata/$line\_1.fastq.gz -o 1-fastp/$line\_clean_1.fastq.gz -I 0-rawdata/$line\_2.fastq.gz -O 1-fastp/$line\_clean_2.fastq.gz -j 1-fastp/$line.json -h 1-fastp/$line.html -z 9 -f 0 -t 0 -T 0 -F 0 -l 75 -p --thread=10
done
##這里我添加了-p參數(shù)嵌器,查看過度表達序列的情況庇谆,也可不加,開啟會增加運行時間

輸出文件如下,兩個版本的報告和輸出結果文件。

Fig.3

報告詳解

報告全部包括這些內容饶囚,點擊每一個條目可以收縮或展開


Fig.4

Summary

Summary中包含了本次分析的情況概要。

Fig.5
  • General

版本號梨与、序列循環(huán)數(shù)瞄崇、質控之前的平均長度、質控之后的平均長度、插入片段的峰值

  • Before filtering

數(shù)據(jù)質控之前的(反應測序質量):總的reads長度大渤、總堿基長度、Q20合格率、Q30合格率、GC含量

  • After filtering

質控之后的:內容同上

  • Filtering result

reads的通過率瞳抓、低質量的reads横蜒、含太多N值的reads

Adapters

這里列出了read1和read2從1到幾十位的adapters的發(fā)生次數(shù),以及其他未列出的接頭數(shù)

Fig.6

Insert size estimation

配對末端重疊分析呆馁,不同長度的Insert在reads中占的比例纺腊,相當于是DNA被打斷后的長度分布次氨。當插入片段大小<30或> 270犀呼,或包含太多錯誤囊骤,則不能被讀取,比如我這里就有33.870706%的不可讀reads)

Fig.7

Before filtering

質控之前的數(shù)據(jù)質量枷踏、堿基含量以及kmer分析等下梢,可直接在網(wǎng)頁上用鼠標拖動放大縮小以及查看具體數(shù)據(jù)細節(jié)岗屏,或進行圖片保存等操作

reads質量

在不同位置上的堿基質量分布暇屋,一般來講質量應 >30 且波動較小為不錯的數(shù)據(jù)


Fig.8

堿基質量

read各個位置上堿基比例分布联予,這個是為了分析堿基的分離程度。何為堿基分離?已知AT配對裳朋,CG配對暖眼,假如測序過程是比較隨機的話(隨機意味著好),那么在每個位置上A和T比例應該差不多丛肢,C和G的比例也應該差不多,如上圖所示猪勇,兩者之間即使有偏差也不應該太大掀泳,最好平均在1%以內马僻,如果過高诗力,除非有合理的原因俺夕,比如某些特定的捕獲測序所致件甥,否則都需要注意是不是測序過程有什么偏差禁灼。

Fig.9

KMER計數(shù)

fastp對5個堿基長度的所有組合的出現(xiàn)次數(shù)進行了統(tǒng)計,然后把它放在了一張表格中斋荞,表格的每一個元素為深背景白字俺驶,背景越深,則表示重復次數(shù)越多塘幅。這樣,一眼望去猪叙,就可以發(fā)現(xiàn)有哪些異常的信息。鼠標可停留在某一具體組合上看出現(xiàn)次數(shù)和平均占比带膀。


Fig.10

過表達序列

提供了這些overrepresented sequence的序列個數(shù)和占比志珍,還提供了他們在測序cycles中的分布情況,這有利于分析各種問題

Fig.11

參考資料:

生信學習筆記:fastp質控處理生成的report結果解讀

fastp: 一款超快速全功能的FASTQ文件自動化質控+過濾+校正+預處理軟件

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
  • 序言:七十年代末垛叨,一起剝皮案震驚了整個濱河市碴裙,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌点额,老刑警劉巖舔株,帶你破解...
    沈念sama閱讀 218,682評論 6 507
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異还棱,居然都是意外死亡载慈,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 93,277評論 3 395
  • 文/潘曉璐 我一進店門珍手,熙熙樓的掌柜王于貴愁眉苦臉地迎上來办铡,“玉大人辞做,你說我怎么就攤上這事」丫撸” “怎么了秤茅?”我有些...
    開封第一講書人閱讀 165,083評論 0 355
  • 文/不壞的土叔 我叫張陵,是天一觀的道長童叠。 經(jīng)常有香客問我框喳,道長,這世上最難降的妖魔是什么厦坛? 我笑而不...
    開封第一講書人閱讀 58,763評論 1 295
  • 正文 為了忘掉前任五垮,我火速辦了婚禮,結果婚禮上杜秸,老公的妹妹穿的比我還像新娘放仗。我一直安慰自己,他們只是感情好撬碟,可當我...
    茶點故事閱讀 67,785評論 6 392
  • 文/花漫 我一把揭開白布诞挨。 她就那樣靜靜地躺著,像睡著了一般呢蛤。 火紅的嫁衣襯著肌膚如雪惶傻。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 51,624評論 1 305
  • 那天顾稀,我揣著相機與錄音达罗,去河邊找鬼。 笑死静秆,一個胖子當著我的面吹牛粮揉,可吹牛的內容都是我干的。 我是一名探鬼主播抚笔,決...
    沈念sama閱讀 40,358評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼扶认,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了殊橙?” 一聲冷哼從身側響起辐宾,我...
    開封第一講書人閱讀 39,261評論 0 276
  • 序言:老撾萬榮一對情侶失蹤,失蹤者是張志新(化名)和其女友劉穎膨蛮,沒想到半個月后叠纹,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 45,722評論 1 315
  • 正文 獨居荒郊野嶺守林人離奇死亡敞葛,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內容為張勛視角 年9月15日...
    茶點故事閱讀 37,900評論 3 336
  • 正文 我和宋清朗相戀三年誉察,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片惹谐。...
    茶點故事閱讀 40,030評論 1 350
  • 序言:一個原本活蹦亂跳的男人離奇死亡持偏,死狀恐怖驼卖,靈堂內的尸體忽然破棺而出,到底是詐尸還是另有隱情鸿秆,我是刑警寧澤酌畜,帶...
    沈念sama閱讀 35,737評論 5 346
  • 正文 年R本政府宣布,位于F島的核電站卿叽,受9級特大地震影響桥胞,放射性物質發(fā)生泄漏。R本人自食惡果不足惜附帽,卻給世界環(huán)境...
    茶點故事閱讀 41,360評論 3 330
  • 文/蒙蒙 一埠戳、第九天 我趴在偏房一處隱蔽的房頂上張望井誉。 院中可真熱鬧蕉扮,春花似錦、人聲如沸颗圣。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,941評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽在岂。三九已至奔则,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間蔽午,已是汗流浹背易茬。 一陣腳步聲響...
    開封第一講書人閱讀 33,057評論 1 270
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留及老,地道東北人抽莱。 一個月前我還...
    沈念sama閱讀 48,237評論 3 371
  • 正文 我出身青樓,卻偏偏與公主長得像骄恶,于是被迫代替她去往敵國和親食铐。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當晚...
    茶點故事閱讀 44,976評論 2 355

推薦閱讀更多精彩內容