因為課題就是做轉(zhuǎn)錄組測序的享怀，所以基礎知識有一些了解添瓷，接下來從數(shù)據(jù)處理部分開始進行筆記鳞贷。
數(shù)據(jù)初步分析：
使用fastqc進行質(zhì)量分析惰聂，這是一款Java軟件庶近，支持多線程鼻种。寫這篇文章的時候版本是v0.11.7。
軟件前期準備：
- 下載方式有兩種：
官網(wǎng)下載好用filezilla導入linux服務器
直接在服務器中wget http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.7.zip
- 接著安裝
unzip fastqc_v0.11.7.zip-->cd FastQC-- >chmod755 fastqc
最后這個chmod有必要說一下投队，這個權限管理命令
chmod 用3個數(shù)字來表達對 用戶（文件或目錄的所有者），用戶組（同組用戶）扒披，其他用戶 的權限：
如：chmod 755 fastqc
數(shù)字7是表達同時具有讀碟案，寫辆亏，執(zhí)行權限：(7 = 4 + 2+ 1)
讀取--用數(shù)字4表示扮叨；
寫入--用數(shù)字2表示甫匹；
執(zhí)行--用數(shù)字1表示; 三者皆否：0
-
設置完權限后抢韭，還需要將FastQC文件夾（這里請注意是文件夾刻恭，而非fastqc這個可執(zhí)行程序）導入環(huán)境變量
echo 'export PATH=/YOUR/FASTQC PATH/:$PATH' >> ~/.bashrc再
source ~/.bashrc 檢查軟件是否安裝成功
fastqc --help出現(xiàn)幫助信息就可以使用啦鞍匾！
后期軟件使用：
基本格式 （各種參數(shù)+多個文件～支持多線程）
fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] seqfile1 .. seqfileN
-o --outdir FastQC生成的報告文件的儲存路徑
--extract 使用這個參數(shù)是讓程序不打包【默認會打包成一個壓縮文件】
-t --threads 選擇程序運行的線程數(shù)，每個線程會占用250MB內(nèi)存（一般與文件數(shù)量一致就好）
-q --quiet 安靜運行模式【不選這個選項梁棠，程序會實時報告運行的狀況】
結(jié)果分析：
- 如果你有自己的轉(zhuǎn)錄組或者其他數(shù)據(jù)，可以現(xiàn)在測試了
- 如果沒有男娄，想學一下這個軟件流程以及如何解讀結(jié)果模闲，可以下載公共數(shù)據(jù)（下載兩個雙端測序文件共4個）
# 順便說一下這里用到了curl -O（保留遠程文件的文件名） -L（對于自動跳轉(zhuǎn)的網(wǎng)頁，curl 就會跳轉(zhuǎn)到新的網(wǎng)址）
# 當然用wget也可以翁授，至于二者區(qū)別贮配，可以參考https://www.cnblogs.com/lsdb/p/7171779.html
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1976948_2.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_1.fastq.gz
curl -O -L https://s3-us-west-1.amazonaws.com/dib-training.ucdavis.edu/metagenomics-scripps-2016-10-12/SRR1977249_2.fastq.gz
# 下載完可以檢查一下數(shù)據(jù)完整性，當然不是必須
md5sum *.gz
# 質(zhì)控四個文件圆存，我們可以采用四線程
# 大概用時 real 0m23.344s （如果你也想統(tǒng)計時間，就在命令前加time）
fastqc *.gz -t 4
#將結(jié)果.html用filezilla導出油讯，瀏覽器查看
-
首先看到的是一個總結(jié)報告，
左邊這一欄會告訴你兔综，哪些是提出警告的（??表示）邻奠，那些是fail的（?）image 接下來一部分一部分解析，共10部分
-
Basic Statistics 基本信息
imageEncoding: 測序平臺編號，現(xiàn)在Sanger/ Illumina 1.8以上都是Phred 33編碼
Total sequences: reads數(shù)量（reads就是高通量測序平臺產(chǎn)生的序列標簽碑隆，翻譯為讀段！）
Sequence length: 測序長度
%GC: GC含量： 需要重點關注劫狠，可以幫助區(qū)別物種呐矾，人類細胞42%左右
Per base sequence quility：每個測序read各堿基質(zhì)量【十分重要懦砂！】
橫軸：測序序列的1-251個堿基蜒犯；
縱軸：質(zhì)量得分，score = -10 * log10（error）荞膘，例如錯誤率error為1%罚随，那么算出的score就是20
箱線圖boxplot：對每一個堿基的質(zhì)量的統(tǒng)計。箱子上面的須（up bar）為90%分位數(shù)羽资，下面的須（down bar）為10%分位數(shù)，箱子中的紅線為中位數(shù)即50%分位數(shù)削罩，箱子頂（upside）為75%分位數(shù)，箱子低（downside）為25%分位數(shù)费奸。這個boxplot的意義：一是看數(shù)據(jù)是否具有對稱性弥激；二是看數(shù)據(jù)分布差異，這里主要利用了第二點愿阐。bar的跨度越大微服，說明數(shù)據(jù)越不穩(wěn)定。
藍色的線將各個堿基的質(zhì)量平均值連接起來
解釋一下：圖中藍線的走勢為何先高后低缨历？因為目前采用的邊合成邊測序使用的是化學方法促使鏈由5'向3'延伸以蕴，也就是利用了DNA聚合酶。剛開始測序辛孵，合成反應還不是很穩(wěn)定丛肮，但是酶的質(zhì)量還很好，所以會在高質(zhì)量區(qū)域內(nèi)有一定的波動（這里的1-30bp）魄缚，后來穩(wěn)定了宝与，但是隨著時間的推移，酶的活力逐漸下降冶匹，特異性也變差习劫，所以越往后出錯幾率越大。
【就像一個司機開車嚼隘，一開始小心謹慎诽里，起步慢，開的也慢飞蛹，慢慢提速谤狡。后來越開越帶勁灸眼，但是也越來越困，疲勞駕駛?cè)菀壮鍪?/strong>】一般能用的數(shù)據(jù)都要求至少Q(mào)20豌汇，也就是下四分位（10%分位數(shù)）的質(zhì)量值要大于20幢炸。因此這里的189bp后面的需要切除，才能繼續(xù)分析
二代測序拒贱，最好是達到Q20的堿基要在95%以上（最差不低于90%）宛徊，Q30要求大于85%（最差也不要低于80%）
- Per sequence quility scores：每條序列 質(zhì)量統(tǒng)計
橫軸：質(zhì)量值0-40，也即是Q值
縱軸：每個質(zhì)量值對應的read數(shù)
我們的例子中一條read有251bp逻澳， 那么其中任意一條的251bp的質(zhì)量平均值就是這條read的質(zhì)量值闸天。只要大部分都高于20說明比較正常
- Per base sequence content：read各個位置堿基比例分布
橫軸：各堿基位置；縱軸：堿基百分比
四條線四種顏色代表四種堿基在每個位置的平均含量（一個位置會測很多reads斜做，然后求一個平均）
一般來講苞氮，A=T， C=G瓤逼， 但是剛開始測序儀不穩(wěn)定可能出現(xiàn)波動笼吟，這是正常的。一般不是波動特別大的霸旗，像這里cut掉前5bp就夠了贷帮。另外如果A、T 或 C诱告、G間出現(xiàn)偏差撵枢，只要在1%以內(nèi)都是可以接受的
- Per sequence GC content： 序列平均GC分布
橫軸為平均GC含量； 縱軸為每個GC含量對應的序列數(shù)量
藍線為系統(tǒng)計算得到的理論分布精居；紅線為測量值锄禽，二者越接近越好
-
這里不相符可能有兩個原因：
前面提到了，GC可以作為物種特異性根據(jù)靴姿，這里出現(xiàn)了其他的峰有可能混入了其他物種的DNA沃但；
目前二代測序基本都會有序列偏向性(所說的 bias)，也就是某些特定區(qū)域會被反復測序佛吓，以至于高于正常水平绽慈，變相說明測序過程不夠隨機。這種現(xiàn)象會對以后的變異檢測以及CNV分析造成影響
如果出現(xiàn)怎么辦辈毯？-- 把和我們使用物種GC-content有差異的reads拿出來做blast坝疼，來確認是否為某些雜菌
- Per base N content： N含量分布
N是指儀器不能識別ATCG時給出的結(jié)果，一般不會出現(xiàn)谆沃。但是如果出現(xiàn)并且量還很大钝凶，應該就是測序系統(tǒng)或者試劑的問題
任意位置的N的比例超過5%，報"WARN"；任意位置的N的比例超過20%耕陷，報"FAIL"
- Sequence length distribution： 序列長度統(tǒng)計
理想情況下闲孤，測得的序列長度應該是相等的锉矢。實際上總有些偏差
當reads長度不一致時報"WARN"崖叫；當有長度為0的read時報“FAIL”
這里顯示大部分都落在251bp這個測序長度上巍耗，有少量為250或252bp，但這不影響嗜诀；如果偏差很大就不可信了
- Sequence duplication level：統(tǒng)計序列完全一樣的reads的頻率
橫坐標是duplication的次數(shù)猾警；縱坐標是duplicated reads的數(shù)目（紅線）
解釋下橫坐標為何會有>10, >50等出現(xiàn)：測序的原始數(shù)據(jù)很大，如果每一條reads都統(tǒng)計隆敢，將耗時很久发皿。這里軟件只采用了數(shù)據(jù)的前200,000條reads統(tǒng)計其在全部數(shù)據(jù)中的重復數(shù)目,另外大于75bp的reads只取50bp進行比較。重復數(shù)大于10的reads被合并統(tǒng)計成了>10拂蝎，以此類推...
unique reads總數(shù)（藍線）作為100%穴墅，上圖中可以看出，大概僅有2%的uniqe reads可以觀察到兩次重復温自。也就是說玄货，我們這里的非unique reads占總數(shù)比例僅有2%左右。
-
正常情況下的確悼泌，測序深度越高誉结，越容易產(chǎn)生一定程度的duplication。高程度的duplication level券躁，提示我們可能有bias的存在（如建庫過程中的PCR duplication）。
另外和做的項目也有關掉盅，一般轉(zhuǎn)錄組測序的結(jié)果中duplication level都比較高也拜，60-70%都正常，這是因為轉(zhuǎn)錄組測的是基因的覆蓋深度趾痘，各個基因表達量不同慢哈，如果某個基因覆蓋度較高【tip：覆蓋度是指基因/轉(zhuǎn)錄組測序測到的部分占整個組的比例】，那么測的部分就越多永票，相對應的duplication也會更高卵贱；對于外顯子組測序來說，一般覆蓋度比較一致侣集，這里出現(xiàn)了duplication就不太正常键俱。
當非unique的reads占總數(shù)的比例大于20%時，報"WARN"世分；當非unique的reads占總數(shù)的比例大于50%時编振，報"FAIL“
- Overrepresented sequences：大量重復序列
和第8個duplication計算一樣，也是取前200,000進行統(tǒng)計臭埋，大于75bp只取50bp踪央。
發(fā)現(xiàn)超過總reads數(shù)0.1%的reads時報”WARN“臀玄，當發(fā)現(xiàn)超過總reads數(shù)1%的reads時報”FAIL“
- Adapter content: 接頭含量
表示序列中兩端adapter的情況
軟件內(nèi)置了四種常用的測序接頭序列， fastqc 有一個參數(shù)-a可以自定義接頭序列
此圖中使用的illumina universal adapter并未去除畅蹂，后期再使用trimmomatic/cutadapt去接頭
我們在學習測序原理的時候知道了健无，接頭的作用一個是能夠使得序列結(jié)合到flowcell上，另外一個多樣本測序時利用接頭旁邊的index加以區(qū)分
什么情況下能夠測到接頭呢液斜？ 一般測序read的長度大于插入片段（待測序列）的長度時會發(fā)生累贤。對于WGS建庫測序來講，一般不會發(fā)生旗唁，因為他的待測序列要幾百bp畦浓，而測序也就150bp算高了。但是對于RNA-seq检疫，一般測序序列比較短讶请，尤其是miRNA，只有幾十bp屎媳，這是就會測到read末尾的接頭
- （還有一類這里沒體現(xiàn)）Kmer content： 重復短序列
表示：在序列中某些特征的短序列重復出現(xiàn)的次數(shù)
-
這個圖是轉(zhuǎn)錄組測序的一個文件夺溢，可以看到6-9bp幾種短序列都出現(xiàn)了好多次。出現(xiàn)的原因可能是：
沒有去除軟件內(nèi)置的adapter或者沒有使用-a參數(shù)自定義adapter
序列本身重復度較高烛谊，例如在建庫PCR過程出現(xiàn)序列偏向性bias--> 這在轉(zhuǎn)錄組測序中確實存在
數(shù)據(jù)過濾：
主要針對接頭序列和低質(zhì)量序列
工具： 有許多工具能干這事风响，例如SOAPnuke、cutadapt丹禀、untrimmed状勤、sickle和seqtk等，其中經(jīng)常用到的是Trimmomatic（也是一個java程序）
-
Trimmomatic:
-
安裝：官網(wǎng)下載我使用的是0.36版本
image下載后直接解壓安裝双泪，其中有一個trimmomatic-0.36.jar是執(zhí)行文件持搜，使用時輸入java -jar trimmomatic-0.36.jar就可以。另外還有一個adapter文件夾焙矛，里面存放了常用的illumina測序儀接頭序列fasta格式,后續(xù)處理接頭需要制定具體文件葫盼。
image
-
-
關于adapter： 目前絕大部分的illumina的Hiseq和Miseq系列使用的都是Truseq3，過去的GA2測序儀使用的是Truseq2村斟，PE/SE對應單端還是雙端測序 贫导， 如果使用的不是illumina測的，可以按照里面的格式自己新建一個接頭文件蟆盹，但其中的命名要注意孩灯。詳情見官網(wǎng)主頁
image 以PE測序為例說一下命令參數(shù)設置：
java -jar <path to trimmomatic.jar> PE [-threads <threads 線程數(shù)]
[-phred33 | -phred64] 質(zhì)量體系，默認64逾滥，但我們現(xiàn)在一般都是33
[-trimlog <logFile>] log文件
<input 1> <input 2> 雙端測序原始fq文件
<paired output 1> <unpaired output 1> 雙端1輸出文件 钱反、 過濾掉的文件
<paired output 2> <unpaired output 2> 雙端2同理
<step >
<詳細講一下step：>
1\. ILLUMINACLIP： 根據(jù)上面qc部分的測試數(shù)據(jù)，我們設置TruSeq2-PE.fa:2:40:15
TruSeq2-PE.fa是接頭序列；
2是比對時接頭序列時所允許的最大錯誤數(shù)面哥；
40是要求PE的兩條read同時和adapter序列比對哎壳，匹配度加起來超40%，那么就認為這對PE reads含 有adapter尚卫，并在對應的位置需要進行切除归榕；
【那么為何不是匹配100%？因為測序時并不是把 adapter全測了吱涉，只測了一部分】
15 指的是只要某條read的某部分與adapter超過了15%的相似度就認為包含刹泄，就要去除
2\. SLIDINGWINDOW： 滑動窗口長度 我們設置為4:20，代表窗口長度為4怎爵，窗口中的平均質(zhì)量值至少為20特石，否 則會開始切除
3\. LEADING，指的是read開頭的堿基是否要被切除的質(zhì)量閾值鳖链，這里設為2
4\. TRAILING姆蘸，指的是read末尾的堿基是否要被切除的質(zhì)量閾值，這里設為2
5\. MINLEN芙委，指的是read被切除后至少需要保留的長度逞敷，如果低于該長度，會被丟掉灌侣，這里設置25
#一個??范例(這個測試數(shù)據(jù)是phred 64的推捐，所以使用默認值就好)：
java -jar trimmomatic-0.36.jar PE -threads 8 \
-trimlog logfile \
reads_1.fq.gz reads_2.fq.gz \
out.read_1.fq.gz out.trim.read_1.fq.gz \
out.read_2.fq.gz out.trim.read_2.fq.gz \ ILLUMINACLIP:/path/Trimmomatic/adapters/TruSeq2-PE.fa:2:40:15 \
SLIDINGWINDOW:4:20 \
LEADING:2 TRAILING:2 \
MINLEN:25</pre>
# 當然，這只是對一個樣本的雙端測序文件進行操作侧啼，那么問題來了：
# 如果你的樣本比較多呢牛柒？還要手動一個個輸入嗎？
# 雖然說vim中使用快速替換 :%s/AA/BB/g 可以全局替換AA成BB痊乾，但還是有點麻煩
# 能不能讓程序來一個自動化呢皮壁？是可以的！可以看作是上面??腳本的改進版～
# 在腳本中構(gòu)建for循環(huán)
# vi trim.sh
#開頭這樣寫而不寫*.gz的目的是避免了樣本名稱的重復
for filename in *_1.fastq.gz
do
#對于filename（%）向右匹配_的全部內(nèi)容(*),然后這部分去掉符喝，留下這之前的
base=${filename%_*}
java -jar trimmomatic-0.36.jar PE \ # 加上反斜線能讓程序整體更清晰
-threads 8 \
-trimlog logfile \
${base}_1.fastq.gz ${base}_2.fastq.gz \
out.${base}_1.fastq.gz out.trim.${base}_1.fastq.gz \
out.${base}_2.fastq.gz out.trim.${base}_2.fastq.gz \
ILLUMINACLIP:/home/u1239/biosoft/Trimmomatic-0.36/adapters/TruSeq2-PE.fa:2:30:10 \
SLIDINGWINDOW:5:20 \
LEADING:5 TRAILING:5 \
MINLEN:50
done
-
如果報錯，就說明參數(shù)錯誤例如我在運行的時候就有的問題百思不得解甜孤，敲了好幾遍代碼應該沒錯协饲，但就是報錯，換個服務器就好了缴川，難道是平臺問題茉稠？其實并不是,就是因為一點點小地方,有時只是一個代表換行的反斜線
\image
多樣本質(zhì)控：
fastqc對于一兩個樣品還能吃得消，但樣本多了把夸，我們?nèi)绾瓮瑫r查看對比他們的信息呢而线？
這里就可以使用界面更加優(yōu)美的multiqc軟件了, 他就是fastqc結(jié)果的整合我會從一個初學者角度來一步步進行，其中包括了中間犯的錯誤及改正～，如果你也在運行的過程發(fā)現(xiàn)了類似的錯誤膀篮，可以參考一下嘹狞。
分界線1: 嘗試自己安裝
下載地址：
解壓縮后，會發(fā)現(xiàn)和以前安裝的源代碼文件不同誓竿，他沒有直接顯示可執(zhí)行文件磅网，也沒有配置文件。這是因為multiqc是一個python軟件筷屡，這里先看一下setup.py：
直接使用python setup.py會報錯涧偷，因為可能你的服務器并沒有setuptools，這是python依賴的第三方庫毙死。
先安裝setuptools：
下載地址：
解壓縮完setuptools會看到有這些文件燎潮，這也是一般的python軟件常見的：包括了setup.py，easy_install.py等
# 主要使用setup.py
python setup.py build # 編譯
python setup.py install #安裝
在第二步安裝的時候報錯,原因是不能對/usr/local/lib下的python進行操作扼倘，因為不管/usr/bin還是/usr/local/bin确封，都是可讀不可改，如果自己家目錄環(huán)境變量中查不到python的路徑唉锌，那么安裝過程中會自動調(diào)取更上層目錄的python使用隅肥。這往往會引發(fā)一系列問題。如果要自己更改的話袄简，需要自己home目錄下有python
好吧腥放，那么接下來我們再安裝自己的python，畢竟自己的軟件用起來方便：
# 下載地址
wget https://www.python.org/ftp/python/3.6.5/Python-3.6.5.tgz
# .tgz 就等同tar.gz
# 解壓后會看到配置文件configure绿语，設置成自己的軟件環(huán)境變量即可
./configure --prefix=/YOUR PATH/
make
make install
# 可能結(jié)尾會有報錯秃症，但是不影響使用，出現(xiàn)了可執(zhí)行文件python就成功了
# 將python復制到你的環(huán)境變量中就能直接調(diào)取
# 最后使用which python檢測是否安在了自己的目錄下
自己的python安裝完成了吕粹，進入到解壓好的multiqc文件夾下, python setup.py build 編譯會運行成功种柑，然后python setup.py install
接著，multiqc就會出現(xiàn)在了自己的軟件環(huán)境變量中了匹耕，輸入multiqc就會看到
分界線2: 同一件事聚请，換個玩法
目的：安裝multiqc
途徑：conda自動安裝
緣由：好久不用conda安裝軟件，一直堅持源代碼稳其，因為好掌控驶赏。但是有的軟件依賴的包很多，自己又很難發(fā)現(xiàn)這些潛在的關系既鞠，所以想重新試試conda煤傍，但并不是傻瓜式的使用，而是讓conda聽我的話～授之以魚嘱蛋，不如授之以漁
Here we go蚯姆！
安裝： 安裝過程需要注意
可以自定義新文件夾五续，默認是在home/miniconda3/
??不要將conda添加到PATH中，我們只需要把它視作一個保姆babysitter龄恋，而不要當一個管家housekeeper
安裝完conda疙驾，把miniconda/bin下的conda復制到你的/biotools/PATH（這個環(huán)境變量因人而異）
添加源： 清華源和中科大源都不錯，別忘了再添加一個bioconda源
新建conda專屬下載目錄： 你可以在你的biotools目錄下新建一個conda軟件存放目錄篙挽，例如conda_install荆萤。然后把這個文件夾添加到環(huán)境變量。以后你用conda安裝的所有軟件都存放在這里铣卡，
??和你自己安裝的軟件要區(qū)分開链韭， 然后利用conda install -p /PATH/conda_install multiqc 就可以實現(xiàn)multiqc的安裝了
作者：劉小澤
鏈接：http://www.reibang.com/p/101c14c3a1d2
來源：簡書
著作權歸作者所有。商業(yè)轉(zhuǎn)載請聯(lián)系作者獲得授權煮落，非商業(yè)轉(zhuǎn)載請注明出處敞峭。
