質(zhì)量控制

測序質(zhì)量檢查-FastQC

Fastqc
Fastqc website (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/))

質(zhì)量控制的測序質(zhì)量檢測是通過FastQC軟件實現(xiàn)。fastqc可以不設(shè)置任何參數(shù)運行，這樣會直接在當(dāng)前目錄下生成一個質(zhì)量報告的壓縮文件和文件夾，報告是網(wǎng)頁格式。也可以設(shè)置輸出目錄和是否解壓縮(--noextract)，默認(rèn)設(shè)置會解壓縮。命令如下：

fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN

其中--noextract命令是不解壓縮輸出文件。-t參數(shù)是指定使用線程數(shù)，fastqc似乎并不是并行運算，而是通過線程數(shù)同時執(zhí)行多個程序，比如線程數(shù)指定為4，并不是用4個進(jìn)程去跑一個文件，而是同時跑4個文件，不過4個線程速度提高很大，個人測試感覺10倍速度于2個線程。-q為屏蔽進(jìn)程信息并只輸出錯誤信息，-f參數(shù)為指定輸入文件格式(有bam, sam, fastq可選)

fastqc的結(jié)果在v0.11.5版下共有12項。

1. Basic Statistics
   包含文件名(Filename)、文件類型、總序列數(shù)(Total Sequences)、序列長度(Sequence length)這些基本信息。
2. Per base sequence quality
   序列每個堿基的平均質(zhì)量，越高越好，需要注意會有部分序列在開頭部分質(zhì)量差，所以根據(jù)這個圖在做質(zhì)控時選擇兩端都去低質(zhì)量或只去3'末端。
3. Per tile sequence quality
   新版本增加的功能，不太清楚是干啥的。
4. Per sequence quality scores
   序列平均得分的數(shù)量。右側(cè)越高越好，也就是大多數(shù)序列質(zhì)量都得分在30以上。均值低于27(也就是錯誤率大于0.002)記為警告，均值低于20(錯誤率0.01)記為不合格。
5. Per base sequence content
   顯示堿基比例。正常情況下堿基比例應(yīng)該差不多。AT與CG差異大于10%記為警告，大于20%記為不合格。
6. Per sequence GC content
   每條序列GC含量百分比與模式化的正態(tài)分布GC含量相比較，超過15%記為警告，超過30%為不合格。
7. Per base N content
   每個堿基位置N（未測到或不確定堿基）的比例
8. Sequence Length Distribution
   各序列長度占比
9. Sequence Duplication Levels
   重要序列重復(fù)級別
   理想情況下所有序列應(yīng)該是被隨機打斷了，測序后理論上不太會出現(xiàn)完全相同的序列。但由于PCR擴(kuò)增或者污染等原因會造成一些重復(fù)序列，這里顯示重復(fù)序列的比例。為了節(jié)省內(nèi)存和時間，fastqc僅抽取了前100,000條序列，并只要超過75bp的序列就會被截斷到50bp來分析。fastqc的說明文檔對此進(jìn)行了說明，結(jié)果不影響整體結(jié)果的反應(yīng)。
   所提供的圖像有兩條線來反應(yīng)不同重復(fù)級別，藍(lán)線表示整個序列集合中重復(fù)序列的分布，紅線表示去除重復(fù)序列后的結(jié)果。這是新版本提供的功能，v0.10版本只有重復(fù)序列的程度。
   一個好的結(jié)果應(yīng)該是紅線藍(lán)線最左側(cè)越大越好。通常會在紅色線中看到一個高重復(fù)的峰，但同時在藍(lán)色線上消失，這表明重復(fù)序列并不顯著。如果在藍(lán)色的線中有峰值，說明在大量不同的高度重復(fù)序列，這可能是污染或嚴(yán)重的技術(shù)重復(fù)。
   如果非唯一序列數(shù)超過總數(shù)的50%就會被軟件判斷為不合格。
10. Overrepresented sequence
    過表達(dá)序列。一般認(rèn)為打斷的序列只有很少部分會重復(fù)，如果一段序列重復(fù)達(dá)到總值的0.1%記為警告，超過1%記為不合格。

Trimming and Quality Filtering

根據(jù)結(jié)果去接頭（adapter）、引物（Primary）尾巴（Poly-A）等。必須要去的是接頭。常用的軟件有cutadapt、trim_galore等等。一般用cutadapt，很多去接頭軟件的底層其實也是調(diào)用cutadapt。

Cutadapt

眼科中心服務(wù)器cutadapt 1.9.1版本安裝在c0，c10節(jié)點上，需要提交到這兩個節(jié)點才可以運行，否則很多節(jié)點用的是1.4.1，老版本的問題是功能有限，尤其是對于雙端數(shù)據(jù)不支持（如-A參數(shù)）。cutadapt官網(wǎng)對于Illumina接頭去除的說明如下：

If you have reads containing Illumina TruSeq adapters, follow these steps.

Single-end reads as well as the first reads of paired-end data need to be trimmed with A + the “TruSeq Indexed Adapter”. Use only the prefix of the adapter sequence that is common to all Indexed Adapter sequences:

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -o trimmed.fastq.gz reads.fastq.gz

If you have paired-end data, trim also read 2 with the reverse complement of the “TruSeq Universal Adapter”. The full command-line looks as follows:

cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC  -o trimmed.1.fastq.gz -p trimmed.2.fastq.gz  reads.1.fastq.gz reads.2.fastq.gz

The adapter sequences can be found in the document Illumina TruSeq Adapters De-Mystified.

因此單端數(shù)據(jù)只需要用-a參數(shù)去掉“AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC”就可以了。

按照推薦我雙端數(shù)據(jù)(Pair-End)的命令如下：

#PBS -N cutadapt-HBV 
#PBS -j oe 
#PBS -l nodes=c0:ppn=1 
#PBS -l walltime=5000:00:00 
#PBS -q low  

# set up some parameter

outputdir="/public/users/chentingting/Zoubin/HBV/QC"  

cd /public/users/chentingting/Zoubin/HBV  

cutadapt -a AGATCGGAAGAGC -A AGATCGGAAGAGC -q 30 -m 20 --trim-n -O 10 -o $outputdir/Sample_capsule-1.R1.trimmed.fastq.gz -p $outputdir/Sample_capsule-1.R2.trimmed.fastq.gz Sample_capsule-1.R1.fastq.gz Sample_capsule-1.R2.fastq.gz

其中的參數(shù)說明:
-a 序列 正向接頭序列，單端測序只用這個。
-A 序列 反向接頭序列，雙端情況下設(shè)置。
-q 數(shù)字 表示最低質(zhì)量值，在去接頭前先將低于此數(shù)值的bases去除。如果只設(shè)置一個數(shù)值則從3'末端去除，如果用逗號分割兩個數(shù)值則先去5'末端后去3'末端。一般設(shè)為30。

-q [5'CUTOFF,]3'CUTOFF, --quality-cutoff=[5'CUTOFF,]3'CUTOFF

Trim low-quality bases from 5' and/or 3' ends of each read before adapter removal. Applied to both reads if data is paired. If one value is given, only the 3' end is trimmed. If two comma-separated cutoffs are given, the 5' end is trimmed with the first cutoff, the 3' end with the second.

-m 數(shù)字 表示trim后最短bp低于此數(shù)的reads被拋棄，一般設(shè)為20。

-M 數(shù)字 表示長于此數(shù)字的reads被拋棄，默認(rèn)值不限制。

--max-n=COUNT 拋棄有太多N的reads。COUNT如果設(shè)置為整數(shù)，就是按N的絕對個數(shù)來處理；如果設(shè)置為小數(shù)（0到1之間），就按每條reads中N的百分比來處理。

-O 數(shù)字 表示adapt和序列比對最少overlap的值，高于此值就認(rèn)為是接頭并修剪，默認(rèn)是3，個人設(shè)置至少到5。

-o 目錄 Read1的輸出路徑

-p 目錄 Read2的輸出路徑

根據(jù)fastqc的報告，如果是RNA數(shù)據(jù)尾巴較多的情況，最好再去一次PolyA尾巴，少就不用了。

Trim Galore!

Trim Galore 合并了FastQC和Cutadapt到一個程序中。它的優(yōu)勢在于它可以根據(jù)FastQC分析的個體質(zhì)量對每個reads進(jìn)行修剪。同時可以設(shè)置程序?qū)羟泻蟮男蛄杏肍astQC生成一個統(tǒng)計信息。對雙端序列支持也很好。

選項

• --phred33 使用ASCII+33質(zhì)量得分作為Phred得分標(biāo)準(zhǔn)。默認(rèn)選項

• --fastqc 當(dāng)剪切結(jié)束后用默認(rèn)選項對結(jié)果文件進(jìn)行fastqc分析

• --fastqc_args "<ARGS>" 對fastQC設(shè)置參數(shù)。

• --paired 設(shè)置雙端序列

• --dont_gzip 輸出文件不壓縮

• --gzip 壓縮輸出文件，如果輸入文件是壓縮文件則自動壓縮。

• --length <INT> 設(shè)置低于數(shù)值長度的reads就拋棄掉，默認(rèn)值20.

• -q/--quality <INT> 切除質(zhì)量得分低于設(shè)置值的序列。默認(rèn)值20.

• -a/--adapter <STRING> -a參數(shù)后為切除的接頭序列

• -a2/--adapter2 <STRING> 雙端序列中read2所切除的接頭序列，需配合'--paired'參數(shù)。

• --rrbs 這是Trim Galore最獨特的功能，也是我一開始找到使用這個軟件的原因：特異性處理MspI所處理的RRBS樣本數(shù)據(jù)(識別位點：CCGG)

示例：

trim_galore --phred33 --fastqc --fastqc_args "--noextract"  --paired --retain_unpaired    --dont_gzip  Sample_keratopathy-31R1.fastq.gz Sample_keratopathy-31R2.fastq.gz