【RSS】[3] Filtering FASTQ files & cutadapt用法

在分析了QA報告之后擒权，人們可能想放棄some of the reads based on several criteria;
我們將使用兩個不同的工具來執(zhí)行這些過濾步驟：PRINSEQ 和 fastq mcf;

PRINSEQ提供了廣泛的過濾選項，我們可以在manual中了解它們：

PRINSEQ lite-h

根據(jù)我們從質(zhì)量保證報告中了解到的情況阁谆，我們可以決定應用following filters：

zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt1 \
    -out_bad null \
    -trim_qual_right 30 -min_len 32 \
    -ns_max_p 5 \
    -lc_method dust -lc_threshold 10
    
zcat SRR031714_2.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_2_filt1 \
    -out_bad null \
    -trim_qual_right 30 -min_len 32 \
    -ns_max_p 5 \
    -lc_method dust -lc_threshold 10

練習：我們使用哪些標準來丟棄reads碳抄？
解決方案：用prinseq-lite -h查看參數(shù)，指定標準场绿，去除不符合標準參數(shù)的fastq剖效；

練習：如果我們的文件是phred 64格式的，我們應該使用哪個option裳凸？
答：我們應該使用-phred64選項。正確選擇此選項非常重要劝贸，否則軟件將無法正確解釋質(zhì)量字符串姨谷，這將影響我們的輸出。

練習：過濾完fastq文件后映九，再次fastQC QA報告來決定我們是否對結(jié)果滿意梦湘，你認為我們能解決最初發(fā)現(xiàn)的主要問題嗎？
答案：我們可以使用以下命令在篩選的fastq文件上運行FastQC：

fastqc SRR031714_1_filt1.fastq SRR031714_2_filt1.fastq

然后我們可以打開并檢查報告件甥。

練習：通常在過濾步驟之前和之后捌议，跟蹤fastq文件中可用的讀取次數(shù)也是有用的。你能從FastQC報告中收集這些信息嗎引有？現(xiàn)在想象一下我們正在處理更多的文件瓣颅；你能想出一個更自動化的方法來檢查嗎？
答案：我們可以從FastQC報告的第一個表中獲取每個fastq文件的讀取次數(shù)譬正，如果我們希望自動化流程宫补，還可以執(zhí)行以下命令：

grep '^@SRR' SRR031714_1_filt1.fastq | wc -l
# find the lines starting with "@SRR" and count how many there are
# 5,150,415
# 利用wc指令我們可以計算文件的Byte數(shù)檬姥、字數(shù)、或是列數(shù)
# -l或--lines 只顯示行數(shù)

請注意粉怕，在某些篩選步驟中過于嚴格可能會導致信息的嚴重丟失健民。從一開始就丟棄太多的讀取通常不是一個好的選擇，而且還可以依賴mapping step來丟棄低質(zhì)量的讀取贫贝。

如果使用paired-end data秉犹，篩選步驟可能會變得復雜，因為必須確保兩個fastq文件（每個讀取對一個）以相同的順序包含相同數(shù)量的讀取稚晚。有一些工具可以簡化此任務崇堵，例如fastq mcf。我們只需輸入工具名稱即可打印可用選項列表：

fastq-mcf

我們注意到fastq mcf的主要功能是從fastq文件中刪除adapter蜈彼。因此筑辨，我們需要提供一個帶有adapter序列的fasta文件。在我們的案例中幸逆，我們決定只檢查標準Illumina paired-end adapter：

cat adapters.fa

現(xiàn)在我們已經(jīng)很好地理解了所需的輸入棍辕，我們可以繼續(xù)執(zhí)行該工具，嘗試匹配我們之前使用的PRINSEQ：

fastq-mcf adapters.fa SRR031714_1.fastq.gz SRR031714_2.fastq.gz \
    -o SRR031714_1_filt2.fastq -o SRR031714_2_filt2.fastq \
    -q 30 -P 33 -l 32 --max-ns 1

根據(jù)所使用的filtering options还绘，我們可能最終得到一組具有不同長度的reads楚昭，導致我們可能會在下游分析中遇到一些困難；

因此拍顷，為了方便和safe side抚太，我們可以使用針對所有reads的filtering options；

zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt3 \
    -out_bad null \
    -trim_right 5
    
zcat SRR031714_1.fastq.gz | prinseq-lite \
    -fastq stdin \
    -out_good SRR031714_1_filt3 \
    -out_bad null \
    -trim_right 5

# zcat(選項)(參數(shù))
# zcat命令 用于不真正解壓縮文件昔案，就能顯示壓縮包中文件的內(nèi)容的場合尿贫。

總之，您可以應用許多filtering combinations踏揣，具體選項在很大程度上取決于type of data和后續(xù)分析庆亡。我們建議查看PRINSEQ manual以獲得關(guān)于該主題的詳細概述。

使用cutadapt去除adapter

在生物信息學100個基礎(chǔ)問題 —— 第10題讀懂FastQC報告之a(chǎn)dapter與kmer中捞稿，我們知道又谋，測序結(jié)果中可能會有若干條序列存在adapter的信息，而adapter的信息一般是不在基因組上存在的娱局。所以彰亥，在比對之前如果不把adapter去干凈，我相信你會得到1個非常非常低的mapping rate衰齐。

通常情況下任斋，我們都是使用cutadapt這個軟件進行adapter（接頭）序列的去除。cutadapt這個軟件不但支持單端序列耻涛，還支持雙端序列的切除仁卷，同時還支持gz格式的自動壓縮與解壓縮穴翩。1個常用的切除命令類似：

# 在linux 命令行模式下
cutadapt -a ADAPTER_FWD -A ADAPTER_REV -o out.1.fastq -p out.2.fastq reads.1.fastq reads.2.fastq

# -a是第1個文件的adapter序列
# -A是第2個文件的adapter序列
# -o是第1個輸出文件
# -p是第2個輸出文件
# reads.1.fastq 是第1個輸入文件，也就是雙端測序中的read-1
# reads.2.fastq 是第2個輸入文件锦积，也就是雙端測序中的read-2

cutadapt中-a/-A 參數(shù)與-g/-G參數(shù)分別代表什么意思芒帕？Illumina測序過程中，一般不會用到哪個參數(shù)丰介？
cutadapt可以過濾一些非常短的reads背蟆，請解釋其中-m 參數(shù)是什么意思？為什么要過濾一些非常短的reads哮幢？
在測序的過程中带膀，我們經(jīng)常發(fā)現(xiàn)一些序列的3'端的測序質(zhì)量不太好（如圖2所示），即使去掉adapter以后還是需要把低質(zhì)量的序列再去除1次橙垢，從而保證后續(xù)的mapping質(zhì)量垛叨。cutadapt可以使用一些辦法來去除3'端質(zhì)量不太好的序列。請說明用哪個參數(shù)來設(shè)置相關(guān)的cutoff柜某，并簡要說明cutadapt對read質(zhì)量判斷的策略與方法嗽元。

cutadapt軟件的介紹

這個cutadapt軟件是最常用的去adapter的工具。它是基于Python編寫的一個Python包喂击，還是用bioconda安裝：

conda install -c bioconda cutadapt

關(guān)于adapter的使用匯總

image.png

其實主體還是利用fastqc進行質(zhì)量評估剂癌，然后根據(jù)質(zhì)量評估的信息。進行相應的質(zhì)控翰绊。利用到的軟件是cutadapter 和 fastx-toolkit佩谷，其實還有一些別的軟件可以用，比如說trimmomatic或者是Trim Galore监嗜；cutadapter序列過濾參考

illumina的通用adapter序列是什么谐檀？

Illumina Paired End Adapters(cannot be used for multiplexing)
Top adapter
    5' ACACTCTTTCCCTACACGACGCTCTTCCGATC*T 3'
Bottom adapter
    5' P-GATCGGAAGAGCGGTTCAGCAGGAATGCCGAG 3'
    
**來源**
http://bioinformatics.cvr.ac.uk/blog/illumina-adapter-and-primer-sequences/
不同顏色的圖例代表不同的測序通用的adapter；
如果在當時fastqc分析的時候-a選項沒有內(nèi)容裁奇，則默認使用圖例中的通用adapter序列進行統(tǒng)計桐猬。