質(zhì)量控制
測序質(zhì)量檢查-FastQC
Fastqc
Fastqc website (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/))
質(zhì)量控制的測序質(zhì)量檢測是通過FastQC軟件實(shí)現(xiàn)鉴扫。fastqc可以不設(shè)置任何參數(shù)運(yùn)行拟逮，這樣會(huì)直接在當(dāng)前目錄下生成一個(gè)質(zhì)量報(bào)告的壓縮文件和文件夾，報(bào)告是網(wǎng)頁格式崖堤。也可以設(shè)置輸出目錄和是否解壓縮(--noextract)漓拾，默認(rèn)設(shè)置會(huì)解壓縮听系。命令如下：
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
其中--noextract命令是不解壓縮輸出文件。-t參數(shù)是指定使用線程數(shù)芝此，fastqc似乎并不是并行運(yùn)算，而是通過線程數(shù)同時(shí)執(zhí)行多個(gè)程序因痛，比如線程數(shù)指定為4婚苹，并不是用4個(gè)進(jìn)程去跑一個(gè)文件，而是同時(shí)跑4個(gè)文件鸵膏，不過4個(gè)線程速度提高很大膊升，個(gè)人測試感覺10倍速度于2個(gè)線程。-q為屏蔽進(jìn)程信息并只輸出錯(cuò)誤信息较性，-f參數(shù)為指定輸入文件格式(有bam, sam, fastq可選)
fastqc的結(jié)果在v0.11.5版下共有12項(xiàng)用僧。
- Basic Statistics
包含文件名(Filename)、文件類型赞咙、總序列數(shù)(Total Sequences)责循、序列長度(Sequence length)這些基本信息。 -
Per base sequence quality
序列每個(gè)堿基的平均質(zhì)量攀操，越高越好院仿，需要注意會(huì)有部分序列在開頭部分質(zhì)量差，所以根據(jù)這個(gè)圖在做質(zhì)控時(shí)選擇兩端都去低質(zhì)量或只去3'末端速和。 - Per tile sequence quality
新版本增加的功能歹垫，不太清楚是干啥的。 - Per sequence quality scores
序列平均得分的數(shù)量颠放。右側(cè)越高越好排惨，也就是大多數(shù)序列質(zhì)量都得分在30以上。均值低于27(也就是錯(cuò)誤率大于0.002)記為警告碰凶，均值低于20(錯(cuò)誤率0.01)記為不合格暮芭。 - Per base sequence content
顯示堿基比例鹿驼。正常情況下堿基比例應(yīng)該差不多。AT與CG差異大于10%記為警告辕宏，大于20%記為不合格利职。 - Per sequence GC content
每條序列GC含量百分比與模式化的正態(tài)分布GC含量相比較赏廓，超過15%記為警告，超過30%為不合格。 - Per base N content
每個(gè)堿基位置N（未測到或不確定堿基）的比例 - Sequence Length Distribution
各序列長度占比 - Sequence Duplication Levels
重要序列重復(fù)級(jí)別
理想情況下所有序列應(yīng)該是被隨機(jī)打斷了拥褂，測序后理論上不太會(huì)出現(xiàn)完全相同的序列跳昼。但由于PCR擴(kuò)增或者污染等原因會(huì)造成一些重復(fù)序列碌廓，這里顯示重復(fù)序列的比例机蔗。為了節(jié)省內(nèi)存和時(shí)間，fastqc僅抽取了前100,000條序列魔策，并只要超過75bp的序列就會(huì)被截?cái)嗟?0bp來分析匈子。fastqc的說明文檔對此進(jìn)行了說明，結(jié)果不影響整體結(jié)果的反應(yīng)闯袒。
所提供的圖像有兩條線來反應(yīng)不同重復(fù)級(jí)別虎敦，藍(lán)線表示整個(gè)序列集合中重復(fù)序列的分布，紅線表示去除重復(fù)序列后的結(jié)果政敢。這是新版本提供的功能其徙，v0.10版本只有重復(fù)序列的程度。
一個(gè)好的結(jié)果應(yīng)該是紅線藍(lán)線最左側(cè)越大越好喷户。通常會(huì)在紅色線中看到一個(gè)高重復(fù)的峰唾那，但同時(shí)在藍(lán)色線上消失，這表明重復(fù)序列并不顯著褪尝。如果在藍(lán)色的線中有峰值闹获，說明在大量不同的高度重復(fù)序列，這可能是污染或嚴(yán)重的技術(shù)重復(fù)河哑。
如果非唯一序列數(shù)超過總數(shù)的50%就會(huì)被軟件判斷為不合格避诽。 - Overrepresented sequence
過表達(dá)序列。一般認(rèn)為打斷的序列只有很少部分會(huì)重復(fù)璃谨，如果一段序列重復(fù)達(dá)到總值的0.1%記為警告沙庐，超過1%記為不合格。
Trimming and Quality Filtering
根據(jù)結(jié)果去接頭（adapter）佳吞、引物（Primary）尾巴（Poly-A）等拱雏。必須要去的是接頭。常用的軟件有cutadapt底扳、trim_galore等等铸抑。一般用cutadapt，很多去接頭軟件的底層其實(shí)也是調(diào)用cutadapt衷模。
Cutadapt
眼科中心服務(wù)器cutadapt 1.9.1版本安裝在c0鹊汛，c10節(jié)點(diǎn)上菇爪，需要提交到這兩個(gè)節(jié)點(diǎn)才可以運(yùn)行，否則很多節(jié)點(diǎn)用的是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.gzIf 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.gzThe 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 序列 正向接頭序列，單端測序只用這個(gè)届囚。
-A 序列 反向接頭序列有梆，雙端情況下設(shè)置。
-q 數(shù)字 表示最低質(zhì)量值意系，在去接頭前先將低于此數(shù)值的bases去除泥耀。如果只設(shè)置一個(gè)數(shù)值則從3'末端去除，如果用逗號(hào)分割兩個(gè)數(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的絕對個(gè)數(shù)來處理；如果設(shè)置為小數(shù)（0到1之間）足绅，就按每條reads中N的百分比來處理捷绑。
-O 數(shù)字 表示adapt和序列比對最少overlap的值，高于此值就認(rèn)為是接頭并修剪氢妈，默認(rèn)是3粹污，個(gè)人設(shè)置至少到5。
-o 目錄 Read1的輸出路徑
-p 目錄 Read2的輸出路徑
根據(jù)fastqc的報(bào)告允懂，如果是RNA數(shù)據(jù)尾巴較多的情況厕怜，最好再去一次PolyA尾巴，少就不用了蕾总。
Trim Galore!
Trim Galore 合并了FastQC和Cutadapt到一個(gè)程序中粥航。它的優(yōu)勢在于它可以根據(jù)FastQC分析的個(gè)體質(zhì)量對每個(gè)reads進(jìn)行修剪。同時(shí)可以設(shè)置程序?qū)羟泻蟮男蛄杏肍astQC生成一個(gè)統(tǒng)計(jì)信息生百。對雙端序列支持也很好递雀。
選項(xiàng)
--phred33 使用ASCII+33質(zhì)量得分作為Phred得分標(biāo)準(zhǔn)。默認(rèn)選項(xiàng)
--fastqc 當(dāng)剪切結(jié)束后用默認(rèn)選項(xiàng)對結(jié)果文件進(jìn)行fastqc分析
--fastqc_args "<ARGS>" 對fastQC設(shè)置參數(shù)蚀浆。
--paired 設(shè)置雙端序列
--dont_gzip 輸出文件不壓縮
--gzip 壓縮輸出文件缀程，如果輸入文件是壓縮文件則自動(dòng)壓縮搜吧。
--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最獨(dú)特的功能撩满，也是我一開始找到使用這個(gè)軟件的原因：特異性處理MspI所處理的RRBS樣本數(shù)據(jù)(識(shí)別位點(diǎn)：CCGG)
示例：
trim_galore --phred33 --fastqc --fastqc_args "--noextract" --paired --retain_unpaired --dont_gzip Sample_keratopathy-31R1.fastq.gz Sample_keratopathy-31R2.fastq.gz
