作業(yè)要求

需要用安裝好的sratoolkit把sra文件轉(zhuǎn)換為fastq格式的測(cè)序文件，并且用fastqc軟件測(cè)試測(cè)序文件的質(zhì)量！
作業(yè)，理解測(cè)序reads，GC含量，質(zhì)量值，接頭，index，fastqc的全部報(bào)告，搜索中文教程，并發(fā)在論壇上面。
來(lái)源于生信技能樹(shù)：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

實(shí)驗(yàn)過(guò)程

1.fastq-dump將sra數(shù)據(jù)轉(zhuǎn)換成fastq格式

# 需要將作業(yè)2（http://www.itdecent.cn/p/da377252ee96）中下載的測(cè)序數(shù)據(jù)用工具sratoolkit轉(zhuǎn)換成fastq的格式。
$ for ((i=56;i<=62;i++));do fastq-dump --gzip --split-3 -A ~/disk2/sra/SRR35899$i.sra -O ~/disk2/data/rna-seq;done

fastq-dump 用法：

--gzip 使得輸出的結(jié)果是.gz 的格式
--split-3 對(duì)于PE測(cè)序，輸出的結(jié)果是兩個(gè)_1.fastq.gz
-A| --accession 輸入你的sra文件可以是絕對(duì)路徑，我的數(shù)據(jù)來(lái)源是~disk2/sra/SRR35899$i.sra （ 如果你直接寫(xiě)accession，那么fastq-dump 會(huì)默認(rèn)重新下載數(shù)據(jù)，并且會(huì)放在~/ncbi/public/sra目錄下）
-O 是設(shè)置輸出的目錄

fastq格式：

Fastq格式是一種基于文本的存儲(chǔ)生物序列和對(duì)應(yīng)堿基（或氨基酸）質(zhì)量的文件格式。最初由桑格研究所（Wellcome Trust Sanger Institute）開(kāi)發(fā)出來(lái)，現(xiàn)已成為存儲(chǔ)高通量測(cè)序數(shù)據(jù)的事實(shí)標(biāo)準(zhǔn)。

fastq數(shù)據(jù)格式

如何判斷是Phred64 還是 Phred33 ？
ASCII值小于等于58（相應(yīng)的質(zhì)量得分小于等于25）對(duì)應(yīng)的字符只有在Phred+33的編碼中被使用，所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情況下，ASCII值大于等于74的字符只出現(xiàn)在Phred+64中。如果是最近兩年的測(cè)序數(shù)據(jù)，一般都是Phred33形式的。參考文章：http://blog.csdn.net/huyongfeijoe/article/details/51613827

2.Fastqc 進(jìn)行測(cè)序結(jié)果的質(zhì)控

用法：
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
參數(shù)：
-o 輸出目錄，需自己創(chuàng)建目錄
--(no)extract 是否解壓輸出文件，默認(rèn)是自動(dòng)解壓縮zip文件。加上--noextract不解壓文件。
-f 指定輸入文件的類(lèi)型，支持fastq|bam|sam三種格式的文件，默認(rèn)自動(dòng)識(shí)別。
-t 同時(shí)處理的文件數(shù)目。
-c 是contaminant 文件，會(huì)從中搜索overpresent 序列。

$ mkdir -p ~/disk2/data/QC
$ cd ~/disk2/data/rna-seq
# 將所有的數(shù)據(jù)進(jìn)行質(zhì)控，得到zip的壓縮文件和html文件
$ fastqc -o ~/disk2/data/QC *.fastq.gz

質(zhì)控結(jié)果文件

3.質(zhì)控結(jié)果查看

質(zhì)控結(jié)果有14個(gè)html文件，你可以選擇用瀏覽器打開(kāi)查看最終的QC reports。

• 首先來(lái)大概看一下QC結(jié)果報(bào)告。

  
  
  QC可視化結(jié)果——雙擊html文件，在瀏覽器中直接打開(kāi)

• 左邊是目錄概要，可以點(diǎn)擊想要看的結(jié)果，右邊會(huì)跳轉(zhuǎn)到特定詳細(xì)的可視化結(jié)果。綠色代表“通過(guò)”，黃色代表“警告”，紅色代表“不通過(guò)，失敗”。

  
  
  Summary

• Basic Statistics，基本的數(shù)據(jù)統(tǒng)計(jì)包括文件名，文件類(lèi)型，編碼形式，總的序列數(shù)，質(zhì)量差的序列，序列平均長(zhǎng)度，GC含量。

  
  
  基本數(shù)據(jù)統(tǒng)計(jì)
• Per base sequence quality，每個(gè)read各位置堿基的測(cè)序質(zhì)量。橫軸堿基的位置，縱軸是質(zhì)量分?jǐn)?shù)，Quality score=-10log10p（p代表錯(cuò)誤率），所以當(dāng)質(zhì)量分?jǐn)?shù)為40的時(shí)候，p就是0.0001，質(zhì)量算高了。紅色線代表中位數(shù)，藍(lán)色代表平均數(shù)，黃色是25%-75%區(qū)間，觸須是10%-90%區(qū)間（黃色和觸須我不是特別明白）。若任一位置的下四分位數(shù)低于10或者中位數(shù)低于25，出現(xiàn)“警告”；若任一位置的下四分位數(shù)低于5或者中位數(shù)低于20，出現(xiàn)“失敗，F(xiàn)ail”。
  
  各位置堿基質(zhì)量

• Per tile sequence quality，檢查reads中每一個(gè)堿基位置在不同的測(cè)序小孔之間的偏離度，藍(lán)色代表偏離度小，質(zhì)量好，越紅代表偏離度越大，質(zhì)量越差。

  
  
  偏離度
• Per sequence quality scores，reads質(zhì)量的分布，當(dāng)峰值小于27時(shí)，警告；當(dāng)峰值小于20時(shí)，fail。我的報(bào)告峰值在38。
  
  reads質(zhì)量分布
• Per base sequence content，對(duì)所有reads的每一個(gè)位置，統(tǒng)計(jì)ATCG四種堿基的分布，橫軸為位置，縱軸為堿基含量，正常情況下每個(gè)位置每種堿基出現(xiàn)的概率是相近的，四條線應(yīng)該平行且相近。當(dāng)部分位置堿基的比例出現(xiàn)bias時(shí)，即四條線在某些位置紛亂交織，往往提示我們有overrepresented sequence的污染。本結(jié)果前10個(gè)位置，每種堿基頻率有明顯的差別，說(shuō)明有污染。當(dāng)任一位置的A/T比例與G/C比例相差超過(guò)10%，報(bào)"WARN"；當(dāng)任一位置的A/T比例與G/C比例相差超過(guò)20%，報(bào)"FAIL"。
  
  堿基分布
• Per Sequence GC Content，統(tǒng)計(jì)reads的平均GC含量的分布。紅線是實(shí)際情況，藍(lán)線是理論分布（正態(tài)分布，均值不一定在50%，而是由平均GC含量推斷的）。 曲線形狀的偏差往往是由于文庫(kù)的污染或是部分reads構(gòu)成的子集有偏差（overrepresented reads）。形狀接近正態(tài)但偏離理論分布的情況提示我們可能有系統(tǒng)偏差。偏離理論分布的reads超過(guò)15%時(shí)，報(bào)"WARN"；偏離理論分布的reads超過(guò)30%時(shí)，報(bào)"FAIL"。
  
  reads 平均GC含量分布
• Per base N content，當(dāng)測(cè)序儀器不能辨別某條reads的某個(gè)位置到底是什么堿基時(shí)，就會(huì)產(chǎn)生“N”，統(tǒng)計(jì)N的比率。正常情況下，N值非常小。當(dāng)任意位置的N的比例超過(guò)5%，報(bào)"WARN"；當(dāng)任意位置的N的比例超過(guò)20%，報(bào)"FAIL"。
  
  各位置N的reads比率
• Sequence Length Distribution，reads長(zhǎng)度分布，當(dāng)reads長(zhǎng)度不一致時(shí)報(bào)"WARN"；當(dāng)有長(zhǎng)度為0的read時(shí)報(bào)“FAIL”。
  
  reads 長(zhǎng)度分布
• Sequence Duplication Levels，統(tǒng)計(jì)不同拷貝數(shù)的reads的頻率。測(cè)序深度越高，越容易產(chǎn)生一定程度的duplication，這是正常的現(xiàn)象，但如果duplication的程度很高，就提示我們可能有bias的存在。橫坐標(biāo)是duplication的次數(shù)，縱坐標(biāo)是duplicated reads的數(shù)目，以u(píng)nique reads的總數(shù)作為100%。下圖中，大于10個(gè)重復(fù)的reads占總序列的20%以上，其他依次類(lèi)推。當(dāng)非unique的reads占總數(shù)的比例大于20%時(shí)，報(bào)"WARN"；當(dāng)非unique的reads占總數(shù)的比例大于50%時(shí)，報(bào)"FAIL“。
  
  統(tǒng)計(jì)不同拷貝數(shù)的reads的頻率

• Overrepresented sequences，一條序列的重復(fù)數(shù)，因?yàn)橐粋€(gè)轉(zhuǎn)錄組中有非常多的轉(zhuǎn)錄本，一條序列再怎么多也不太會(huì)占整個(gè)轉(zhuǎn)錄組的一小部分（比如1%），如果出現(xiàn)這種情況，不是這種轉(zhuǎn)錄本巨量表達(dá)，就是樣品被污染。這個(gè)模塊列出來(lái)大于全部轉(zhuǎn)錄組1%的reads序列，但是因?yàn)橛玫氖乔?00,000條，所以其實(shí)參考意義不大，完全可以忽略。

  
  
  一條序列的重復(fù)數(shù)

• Adapter content，接頭含量

  
  
  接頭含量

• Kmer content

  
  
  Kmer含量

參考資料：https://www.plob.org/article/5987.html；http://blog.sina.com.cn/s/blog_1319a10ee0102vfbx.html。

4.質(zhì)控結(jié)果批量查看神器——MultiQC

知乎青山屋主寫(xiě)的為知筆記介紹了multiQC軟件--批量顯示QC結(jié)果。

# 利用Anaconda來(lái)安裝MultiQC非常方便，首先得安裝Anaconda，用清華源下載，特別快，而官網(wǎng)實(shí)則難以接受。
# 清華源地址：https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
# 官網(wǎng)：https://www.continuum.io/downloads/
$ cd ~/src
$ wget  https://mirrors.tuna.tsinghua.edu.cn/anaconda/archive/
# 下載到的是shell腳本文件，直接運(yùn)行，安裝完成
$ bash Anaconda2-4.4.0-Linux-x86_64.sh
# 添加 Anaconda Python 免費(fèi)倉(cāng)庫(kù)
$ conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free/
$ conda config --set show_channel_urls yes
# 然后直接安裝MultiQC
$ conda install -c bioconda multiqc
# 測(cè)試
$ multiqc --help
# 進(jìn)入存放QC結(jié)果的文件夾，并執(zhí)行multiqc
$ cd ~/disk2/data/QC
# 掃描結(jié)果文件，忽略html文件
$ multiqc /data/*fastqc.zip --ignore *.html
# 最后會(huì)默認(rèn)生成一個(gè)名為multiqc_report.html文件，用瀏覽器查看，具體看青山屋主的介紹。

參考資料：https://mirrors.tuna.tsinghua.edu.cn/help/anaconda/
https://www.continuum.io/downloads/
http://fbb84b26.wiz03.com/share/s/3XK4IC0cm4CL22pU-r1HPcQQ1iRTvV2GwkwL2AaxYi2fXHP7