今天將會用到第二天建立的rnaseq環(huán)境及對應(yīng)文件路徑和Day3下載好的.sra數(shù)據(jù)
今天主要學(xué)習(xí)內(nèi)容是配置相關(guān)軟件環(huán)境，開始數(shù)據(jù)轉(zhuǎn)換队丝，質(zhì)量控制，了解參考基因組及注釋文件
數(shù)據(jù)轉(zhuǎn)換：其實(shí)這一步就是將下載好的.sra文件轉(zhuǎn)換成fq文件秧饮，因?yàn)檠芯空咝枰獙⑾聶C(jī)的fq文件壓縮成sra上傳到ncbi野芒，我們這一步的工作就是再將其還原呜师。主要用到的工具是fastq-dump娶桦，它和Day3天的perfetch同樣數(shù)據(jù)sratools工具包贾节，所以直接調(diào)用就可以汁汗。
數(shù)據(jù)質(zhì)控：這一步的目的就是通過查看有無接頭和低質(zhì)量堿基等，得知原始數(shù)據(jù)質(zhì)量咋樣栗涂，不好就要進(jìn)行處理知牌，處理之后再不好那可能就是要重新檢測，主要用到的軟件是fastqc和multiqc斤程，其中當(dāng)我們分析數(shù)據(jù)存在多個(gè)樣本時(shí)角寸，為了方便同時(shí)查看質(zhì)量信息我們將使用multiqc進(jìn)行多樣本合并呈現(xiàn)結(jié)果菩混。
數(shù)據(jù)過濾：通過數(shù)據(jù)質(zhì)控知道了數(shù)據(jù)所存在的問題，那么接下來我們要開始對原始數(shù)據(jù)進(jìn)行過濾扁藕，因?yàn)榻宇^盒低質(zhì)量的堿基都會影響后續(xù)的比對沮峡、定量準(zhǔn)確性，及下游分析亿柑。主要用到的工具為trim_galore或者trimmomatic邢疙，SOAPfilter（這三個(gè)區(qū)別）
比對：比對的軟件有很多，主要有三種類型：基于基因組比對：star望薄，hisat2疟游；基于轉(zhuǎn)錄組比對：bowtie，bwa痕支；不基于比對：salmon
定量：推薦使用featureCounts颁虐，它是subread軟件下的一個(gè)小軟件
開始！
第一步-軟件環(huán)境
依然是使用conda進(jìn)行配置卧须，主要軟件包括：數(shù)據(jù)轉(zhuǎn)換另绩，質(zhì)量控制，數(shù)據(jù)過濾花嘶，比對和定量
#激活之前設(shè)置好的rnaseq環(huán)境
conda activate rnaseq 或 source activate rnaseq
#安裝所需軟件; -y:默認(rèn)安裝過程全部yes板熊，當(dāng)conda把所有軟件都安裝好之后才會給反饋
conda install fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon -y
請注意軟件名稱大小寫，很重要察绷，錯(cuò)了容易報(bào)錯(cuò)干签；同時(shí)還有可能存在軟件依賴性，比如這里的featurecounts找不到并報(bào)錯(cuò)拆撼，此處的意思是找遍了所有conda下載channel也沒找到它容劳，這是我們就網(wǎng)頁搜索一下，根據(jù)正確答案將代碼替換就ok啦
當(dāng)我將fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon這些軟件都安裝好后闸度，每一個(gè)軟件我都檢測一下到底安裝好沒竭贩，通過：輸入軟件名直接回車、軟件名 --help莺禁、軟件名 -h留量、軟件名 -help，以上方法均為查看該軟件幫助文件哟冬，當(dāng)幫助文件顯示后則表示為安裝成功楼熄。
但！請看（學(xué)習(xí)不認(rèn)真就要入深坑）
這兩個(gè)軟件沒有顯示幫助文件浩峡，可能是安裝出了問題可岂，那么我通過網(wǎng)頁搜索其安裝方法，進(jìn)行重新安裝翰灾，如下缕粹，但還是出現(xiàn)問題：
問題一：
我以為的真的不是我以為的稚茅，我以為所有的軟件的幫助文件都是這樣看的，結(jié)果我錯(cuò)了平斩，首先subread不是單獨(dú)的軟件亚享，它還包括subread-align，subjunc绘面，其官網(wǎng)：(http://subread.sourceforge.net/)虹蒋，下這個(gè)軟件的目的是為了使用其下面的featureCounts，所以我們直接輸入 featureCounts回車飒货，也可以證明我們這個(gè)軟件安裝成功
問題二：
trim-galore也不顯示幫助文件魄衅，也重新安裝了，沒報(bào)錯(cuò)塘辅，去安裝的路徑（home/miniconda3/envs/rnaseq/bin）下看它也在里面呀晃虫，那這又咋了啊扣墩！真讓人頭大哲银！
我的天！這是我們可以vi ~/.bashrc呻惕，將其添加到變量中export PATH=~/miniconda3/envs/rnaseq/bin:$PATH荆责，source ~/.bashrc激活一下，只是再trim_galore回車就會顯示幫助文件信息了
所以基于以上亚脆，當(dāng)我們遇到安裝軟件之后不顯示幫助文件信息時(shí)處理方式有三種（目前我遇到的做院，誰知道以后還會遇見啥呢）
1.搜索conda 軟件名：根據(jù)網(wǎng)站提示重新安裝
2.通過1沒有報(bào)錯(cuò)，安裝成功濒持，但還是不顯示幫助文件键耕，那我們就要了解這個(gè)軟件是不是唯一的，其內(nèi)部會不會還有其他軟件（類似問題一）
3.查看軟件安裝路徑內(nèi)是否有該軟件柑营，有的話直接添加變量就OK了（類似問題二）
Ok屈雄！配置相關(guān)軟件環(huán)境完成
第二步-數(shù)據(jù)轉(zhuǎn)換
由于目前測序都是雙末端測序，因此將sra轉(zhuǎn)換成fq后會出現(xiàn)兩個(gè)fq文件官套，read1和read2酒奶，將下列腳本寫入名字為sra2fq.sh中（vi sra2fq.sh）
cat srr.ids | while read I ;do
echo $i
#time fastq-dump –gzip –split-3 -A $i 路徑/${i}.sra -O 輸出路徑 &；
done
首先在做循環(huán)腳本時(shí)奶赔，真正要運(yùn)行的命令先用#注釋掉（小澤提的好建議）惋嚎，用echo $i可以先看看這個(gè)循環(huán)是否是可行的，運(yùn)行一下纺阔，看結(jié)果：
確認(rèn)是我們要分析的數(shù)據(jù)列表瘸彤，好了，這時(shí)回到腳本將注釋打開
當(dāng)我們在使用某一個(gè)軟件時(shí)笛钝，幫助軟件時(shí)非常有用的质况，因?yàn)楹芏鄥?shù)我們沒有必要都要記住，可能也記不住玻靡，所以我們要充分利用好幫助信息结榄，查看某個(gè)軟件的幫助信息的命令前面已經(jīng)說了太多了（-help）
如何可以讓一個(gè)安裝好的軟件被我們所用，也就是將它運(yùn)行起來囤捻，這里用featureCounts舉例（因?yàn)槲医K于弄明白了subread的幫助文件怎么顯示了）臼朗；主要三因素，主程序+選項(xiàng)+參數(shù)：featureCounts [option] -a <annotation_file>蝎土，[]里內(nèi)容表示可選视哑，<>里內(nèi)容為必選選項(xiàng)，不可忽略誊涯，在寫命令時(shí)是不需要寫上[]和<>挡毅。
此外，在幫助文件中你會看到有一個(gè)連字符+大/小寫字母暴构，和兩個(gè)連字符+單詞跪呈，兩個(gè)連字符+單詞要比一個(gè)連字符+大/小寫字母說的更詳細(xì)。這兩種方式主要是對選項(xiàng)的調(diào)用取逾，后面才是選項(xiàng)的詳細(xì)描述內(nèi)容耗绿。
數(shù)據(jù)轉(zhuǎn)換命令詳述：
time fastq-dump –gzip –split-3 -A $i 路徑/${i}.sra -O 輸出路徑 1>sra2fq.log 2>&1 &
#time:表示計(jì)時(shí)，就是轉(zhuǎn)換這個(gè)過程用了多久砾隅，可有可無
#fastq-dump:主程序
#--gzip：表示將解壓出來的fq還是壓縮成gzip格式误阻，目的是節(jié)省空間內(nèi)存
#--split-3：https://www.biostars.org/p/156909/，雖然其中有個(gè)3晴埂，但其實(shí)結(jié)果一般也就出2個(gè)文件
#-A指定輸出結(jié)果名
#i：SRRxxxx
#-O：指出輸出文件路徑
#1>sra2fq.log：將結(jié)果的正確日志文件輸出到sra2fq.log
#2>&1:2表示將錯(cuò)誤信息堕绩，并將其追加到1的正確日志中。也可以將2錯(cuò)誤信息丟棄到linux黑洞中去：2>/dev/null
#&：實(shí)現(xiàn)并行運(yùn)算的意思邑时，說白了就是同時(shí)跑這個(gè)循環(huán)奴紧，對i進(jìn)行同時(shí)數(shù)據(jù)轉(zhuǎn)換
--split-3說白了就是輸出3個(gè)結(jié)果文件如下，但是我們往下分析的時(shí)候只用read1和read2晶丘，所以如果只想要2個(gè)結(jié)果文件黍氮，可以用--split-files；
好了浅浮，現(xiàn)在可以投任務(wù)了 nohup sh sra2fq.sh &沫浆，投任務(wù)這個(gè)時(shí)間可以準(zhǔn)備下載參考基因組和注釋文件了
第三步-參考文件
為什么需要參考文件？都有哪些參考文件滚秩？
第一個(gè)參考文件：參考基因組
轉(zhuǎn)錄組最終的結(jié)果就是要比較表達(dá)量的差異专执，那么一堆隨機(jī)排列的堿基怎么看表達(dá)量差異呢？就需要一個(gè)“標(biāo)準(zhǔn)”郁油，我們測到的不同數(shù)據(jù)都和這個(gè)標(biāo)準(zhǔn)去做比較本股，比上的多攀痊，我們就認(rèn)為表達(dá)量高，相反表達(dá)量少拄显，對于比對不上的苟径，我們認(rèn)為不表達(dá)。其實(shí)這里說的“標(biāo)準(zhǔn)”也就是第一個(gè)參考文件指的就是參考基因組躬审，人類基因組經(jīng)歷了從hg16棘街、hg17、hg18承边、hg19遭殉、hg38版本，最新版本即為hg38博助。
人類參考基因組：gz壓縮文件800多M险污，解壓后3G。
有3種主流的數(shù)據(jù)庫可供使用下載不同物種基因組數(shù)據(jù)：NCBI翔始、UCSC罗心、Ensembl；
-NCBI：ftp://ftp.ncbi.nlm.nih.gov/genomes/
-UCSC：http://hgdownload.cse.ucsc.edu/downloads.html
-Ensembl：ftp://ftp.ensembl.org/pub/
第二個(gè)參考文件：基因組注釋文件
有了基因組城瞎，可以讓我們?nèi)ヌ剿鞯玫降臄?shù)據(jù)reads中哪些序列匹配到基因組上的對應(yīng)位置渤闷，但是我們還是不知道對應(yīng)的是什么，如果我們把基因組比作“尋路”的過程脖镀，那么注釋文件就是給找好的路指定“路牌”飒箭，主要使用兩個(gè)數(shù)據(jù)庫：Gencode：Ensembl，下載得到的文件主要是GTF和GFF：
Gencode：https://www.gencodegenes.org/蜒灰，包括人和小鼠的
Ensembl：ftp://ftp.ensembl.org/pub/current_gtf
關(guān)于GTF和GFF文件的信息請學(xué)習(xí)：小澤分享的GTF與GFF的小趣事
那么現(xiàn)在開始下載兩個(gè)參考文件吧弦蹂！
請將兩個(gè)參考文件下載到之前創(chuàng)建的ref文件目錄下（cd ref）
#genome（從ensembl下載）
wget -c ftp://ftp.ensembl.org/pub/current_fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
#annotation（從ensembl下載）
wget -c ftp://ftp.ensembl.org/pub/current_gtf/homo_sapiens/Homo_sapiens.GRCh38.96.gtf.gz
# -c的含義是斷點(diǎn)續(xù)傳
# 如果下載速度慢，可以放后臺
#解壓
gzip -d Homo_sapiens.GRCh38.96.gtf.gz
gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
第四步-利用fastqc進(jìn)行質(zhì)控
通過fastq-dump將.sra文件轉(zhuǎn)換成fq后强窖，共生成8對SRRxxxxxx_1.fastq.gz和SRRxxxxxx_2.fastq.gz凸椿，接下來對其進(jìn)行fastqc質(zhì)控，在qc文件目錄下進(jìn)行（cd qc）
fastqc 路徑/raw/*.gz -o ./ -t 10
# -t --threads 選擇程序運(yùn)行的線程數(shù)翅溺，越多越快哦脑漫；-o輸出路徑；*.gz表示該目錄下所有的.gz文件咙崎；
因?yàn)槲业膄astq-dump命令還沒有跑完优幸，所以我單獨(dú)跑已完成的.fastq.gz（此處我跑兩個(gè)SRRxxxxx數(shù)據(jù)文件，方便后面multiqc使用及結(jié)果演示）褪猛，如下：
運(yùn)行結(jié)果會按照每個(gè)fq文件以html格式輸出网杆，然后利用filezilla導(dǎo)出結(jié)果（我用的導(dǎo)出結(jié)果是WinSCP）如下：
右鍵.html即可查看結(jié)果，如下：
至于這個(gè)Report怎么解讀，請學(xué)習(xí)繼續(xù)小澤的優(yōu)秀文章====不要"太重視“fastqc的結(jié)果
FastQC Report詳解：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/
因?yàn)槲覀冏罱K會輸出16個(gè).html碳却，一個(gè)一個(gè)看實(shí)在是太麻煩队秩，所以我們利用multiqc可以將結(jié)果合并呈現(xiàn)
額外說明一下multiqc安裝有兩種方法：官網(wǎng)：https://multiqc.info/docs/
# pip安裝
pip install git+https://github.com/ewels/MultiQC.git #Installation with pip
# conda安裝
conda install -c bioconda multiqc # Installing with conda
插播：在等待安裝multiqc時(shí)，所有數(shù)據(jù)fastqc結(jié)束追城，如下刹碾；
運(yùn)行MultiQC
直接指定MultiQC要分析的文件路徑即可燥撞，若數(shù)據(jù)在當(dāng)前目錄下輸入multiqc .即可座柱。
multiqc .
#使用--ignore忽略掉某些文件
#multiqc . --ignore *_R2*
multiqc結(jié)果解析，請學(xué)習(xí)：http://www.reibang.com/p/f83626fd1fa1
坑物舒！坑色洞！關(guān)于multiqc
我以為的又不是我以為的，2個(gè)小時(shí)出坑
一定要清楚multiqc是干啥的冠胯，我一直以為是和fastqc一樣的進(jìn)行質(zhì)控火诸，錯(cuò)！我以為錯(cuò)了荠察，其實(shí)它是將fastqc質(zhì)檢后生成的.zip文件進(jìn)行質(zhì)量信息結(jié)果合并呈現(xiàn)置蜀，所以我們應(yīng)該在qc這個(gè)路徑下對.zip文件進(jìn)行multiqc .
頓悟這個(gè)問題的心里路程：
當(dāng)我第一次運(yùn)行multiqc .最先出現(xiàn)這個(gè)報(bào)錯(cuò)，下面繼續(xù)出錯(cuò)顯示W(wǎng)ARNING
頓時(shí)蒙圈了悉盆，別人也沒這個(gè)情況呀盯荤，我就開始搜索，搜呀搜呀搜呀搜焕盟，第一個(gè)警告解決秋秤，幫助來自https://github.com/yaml/pyyaml/wiki/PyYAML-yaml.load(input)-Deprecation#footnotes
找到config.py文件，對其進(jìn)行編輯脚翘，vi 路徑/config.py
在vi編輯器中查找yaml.load灼卢，一般都是yaml.load（f），這是我們要在其后面加上Loader=yaml.FullLoader来农，把所有的yaml.load中都加上鞋真，之后保存wq，退出vi
yaml.load (f, Loader=yaml.FullLoader)
接下來我們開始解決WARNING沃于，用了1個(gè)小時(shí)解決的涩咖，我都想爆粗口，現(xiàn)在想想就一句話揽涮，multiqc處理的是fastqc質(zhì)控后生成的.zip抠藕，找對路徑找對文件執(zhí)行就完事了，要被自己蠢哭了蒋困，最后生成multiqc_data和.html盾似，使用WinSCP打開.html
