今天將會用到第二天建立的rnaseq環(huán)境及對應文件路徑和Day3下載好的.sra數(shù)據(jù)
今天主要學習內(nèi)容是配置相關(guān)軟件環(huán)境，開始數(shù)據(jù)轉(zhuǎn)換，質(zhì)量控制，了解參考基因組及注釋文件
數(shù)據(jù)轉(zhuǎn)換：其實這一步就是將下載好的.sra文件轉(zhuǎn)換成fq文件，因為研究者需要將下機的fq文件壓縮成sra上傳到ncbi，我們這一步的工作就是再將其還原。主要用到的工具是fastq-dump，它和Day3天的perfetch同樣數(shù)據(jù)sratools工具包，所以直接調(diào)用就可以。
數(shù)據(jù)質(zhì)控：這一步的目的就是通過查看有無接頭和低質(zhì)量堿基等，得知原始數(shù)據(jù)質(zhì)量咋樣，不好就要進行處理，處理之后再不好那可能就是要重新檢測，主要用到的軟件是fastqc和multiqc，其中當我們分析數(shù)據(jù)存在多個樣本時，為了方便同時查看質(zhì)量信息我們將使用multiqc進行多樣本合并呈現(xiàn)結(jié)果。
數(shù)據(jù)過濾：通過數(shù)據(jù)質(zhì)控知道了數(shù)據(jù)所存在的問題，那么接下來我們要開始對原始數(shù)據(jù)進行過濾，因為接頭盒低質(zhì)量的堿基都會影響后續(xù)的比對、定量準確性，及下游分析。主要用到的工具為trim_galore或者trimmomatic，SOAPfilter（這三個區(qū)別）
比對：比對的軟件有很多，主要有三種類型：基于基因組比對：star，hisat2；基于轉(zhuǎn)錄組比對：bowtie，bwa；不基于比對：salmon
定量：推薦使用featureCounts，它是subread軟件下的一個小軟件
開始！
第一步-軟件環(huán)境
依然是使用conda進行配置，主要軟件包括：數(shù)據(jù)轉(zhuǎn)換，質(zhì)量控制，數(shù)據(jù)過濾，比對和定量
#激活之前設置好的rnaseq環(huán)境
conda activate rnaseq 或 source activate rnaseq
#安裝所需軟件; -y:默認安裝過程全部yes，當conda把所有軟件都安裝好之后才會給反饋
conda install fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon -y
請注意軟件名稱大小寫，很重要，錯了容易報錯；同時還有可能存在軟件依賴性，比如這里的featurecounts找不到并報錯，此處的意思是找遍了所有conda下載channel也沒找到它，這是我們就網(wǎng)頁搜索一下，根據(jù)正確答案將代碼替換就ok啦
當我將fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon這些軟件都安裝好后，每一個軟件我都檢測一下到底安裝好沒，通過：輸入軟件名直接回車、軟件名 --help、軟件名 -h、軟件名 -help，以上方法均為查看該軟件幫助文件，當幫助文件顯示后則表示為安裝成功。
但！請看（學習不認真就要入深坑）
這兩個軟件沒有顯示幫助文件，可能是安裝出了問題，那么我通過網(wǎng)頁搜索其安裝方法，進行重新安裝，如下，但還是出現(xiàn)問題：
問題一：
我以為的真的不是我以為的，我以為所有的軟件的幫助文件都是這樣看的，結(jié)果我錯了，首先subread不是單獨的軟件，它還包括subread-align，subjunc，其官網(wǎng)：(http://subread.sourceforge.net/)，下這個軟件的目的是為了使用其下面的featureCounts，所以我們直接輸入 featureCounts回車，也可以證明我們這個軟件安裝成功
問題二：
trim-galore也不顯示幫助文件，也重新安裝了，沒報錯，去安裝的路徑（home/miniconda3/envs/rnaseq/bin）下看它也在里面呀，那這又咋了?。≌孀屓祟^大！
我的天！這是我們可以vi ~/.bashrc，將其添加到變量中export PATH=~/miniconda3/envs/rnaseq/bin:$PATH，source ~/.bashrc激活一下，只是再trim_galore回車就會顯示幫助文件信息了
所以基于以上，當我們遇到安裝軟件之后不顯示幫助文件信息時處理方式有三種（目前我遇到的，誰知道以后還會遇見啥呢）
1.搜索conda 軟件名：根據(jù)網(wǎng)站提示重新安裝
2.通過1沒有報錯，安裝成功，但還是不顯示幫助文件，那我們就要了解這個軟件是不是唯一的，其內(nèi)部會不會還有其他軟件（類似問題一）
3.查看軟件安裝路徑內(nèi)是否有該軟件，有的話直接添加變量就OK了（類似問題二）
Ok！配置相關(guān)軟件環(huán)境完成
第二步-數(shù)據(jù)轉(zhuǎn)換
由于目前測序都是雙末端測序，因此將sra轉(zhuǎn)換成fq后會出現(xiàn)兩個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)腳本時，真正要運行的命令先用#注釋掉（小澤提的好建議），用echo $i可以先看看這個循環(huán)是否是可行的，運行一下，看結(jié)果：
確認是我們要分析的數(shù)據(jù)列表，好了，這時回到腳本將注釋打開
當我們在使用某一個軟件時，幫助軟件時非常有用的，因為很多參數(shù)我們沒有必要都要記住，可能也記不住，所以我們要充分利用好幫助信息，查看某個軟件的幫助信息的命令前面已經(jīng)說了太多了（-help）
如何可以讓一個安裝好的軟件被我們所用，也就是將它運行起來，這里用featureCounts舉例（因為我終于弄明白了subread的幫助文件怎么顯示了）；主要三因素，主程序+選項+參數(shù)：featureCounts [option] -a <annotation_file>，[]里內(nèi)容表示可選，<>里內(nèi)容為必選選項，不可忽略，在寫命令時是不需要寫上[]和<>。
此外，在幫助文件中你會看到有一個連字符+大/小寫字母，和兩個連字符+單詞，兩個連字符+單詞要比一個連字符+大/小寫字母說的更詳細。這兩種方式主要是對選項的調(diào)用，后面才是選項的詳細描述內(nèi)容。
數(shù)據(jù)轉(zhuǎn)換命令詳述：
time fastq-dump –gzip –split-3 -A $i 路徑/${i}.sra -O 輸出路徑 1>sra2fq.log 2>&1 &
#time:表示計時，就是轉(zhuǎn)換這個過程用了多久，可有可無
#fastq-dump:主程序
#--gzip：表示將解壓出來的fq還是壓縮成gzip格式，目的是節(jié)省空間內(nèi)存
#--split-3：https://www.biostars.org/p/156909/，雖然其中有個3，但其實結(jié)果一般也就出2個文件
#-A指定輸出結(jié)果名
#i：SRRxxxx
#-O：指出輸出文件路徑
#1>sra2fq.log：將結(jié)果的正確日志文件輸出到sra2fq.log
#2>&1:2表示將錯誤信息，并將其追加到1的正確日志中。也可以將2錯誤信息丟棄到linux黑洞中去：2>/dev/null
#&：實現(xiàn)并行運算的意思，說白了就是同時跑這個循環(huán)，對i進行同時數(shù)據(jù)轉(zhuǎn)換
--split-3說白了就是輸出3個結(jié)果文件如下，但是我們往下分析的時候只用read1和read2，所以如果只想要2個結(jié)果文件，可以用--split-files；
好了，現(xiàn)在可以投任務了 nohup sh sra2fq.sh &，投任務這個時間可以準備下載參考基因組和注釋文件了
第三步-參考文件
為什么需要參考文件？都有哪些參考文件？
第一個參考文件：參考基因組
轉(zhuǎn)錄組最終的結(jié)果就是要比較表達量的差異，那么一堆隨機排列的堿基怎么看表達量差異呢？就需要一個“標準”，我們測到的不同數(shù)據(jù)都和這個標準去做比較，比上的多，我們就認為表達量高，相反表達量少，對于比對不上的，我們認為不表達。其實這里說的“標準”也就是第一個參考文件指的就是參考基因組，人類基因組經(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/
第二個參考文件：基因組注釋文件
有了基因組，可以讓我們?nèi)ヌ剿鞯玫降臄?shù)據(jù)reads中哪些序列匹配到基因組上的對應位置，但是我們還是不知道對應的是什么，如果我們把基因組比作“尋路”的過程，那么注釋文件就是給找好的路指定“路牌”，主要使用兩個數(shù)據(jù)庫：Gencode：Ensembl，下載得到的文件主要是GTF和GFF：
Gencode：https://www.gencodegenes.org/，包括人和小鼠的
Ensembl：ftp://ftp.ensembl.org/pub/current_gtf
關(guān)于GTF和GFF文件的信息請學習：小澤分享的GTF與GFF的小趣事
那么現(xiàn)在開始下載兩個參考文件吧！
請將兩個參考文件下載到之前創(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的含義是斷點續(xù)傳
# 如果下載速度慢，可以放后臺
#解壓
gzip -d Homo_sapiens.GRCh38.96.gtf.gz
gzip -d Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
第四步-利用fastqc進行質(zhì)控
通過fastq-dump將.sra文件轉(zhuǎn)換成fq后，共生成8對SRRxxxxxx_1.fastq.gz和SRRxxxxxx_2.fastq.gz，接下來對其進行fastqc質(zhì)控，在qc文件目錄下進行（cd qc）
fastqc 路徑/raw/*.gz -o ./ -t 10
# -t --threads 選擇程序運行的線程數(shù)，越多越快哦；-o輸出路徑；*.gz表示該目錄下所有的.gz文件；
因為我的fastq-dump命令還沒有跑完，所以我單獨跑已完成的.fastq.gz（此處我跑兩個SRRxxxxx數(shù)據(jù)文件，方便后面multiqc使用及結(jié)果演示），如下：
運行結(jié)果會按照每個fq文件以html格式輸出，然后利用filezilla導出結(jié)果（我用的導出結(jié)果是WinSCP）如下：
右鍵.html即可查看結(jié)果，如下：
至于這個Report怎么解讀，請學習繼續(xù)小澤的優(yōu)秀文章====不要"太重視“fastqc的結(jié)果
FastQC Report詳解：http://www.bioinformatics.babraham.ac.uk/projects/fastqc/Help/
因為我們最終會輸出16個.html，一個一個看實在是太麻煩，所以我們利用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ù)據(jù)fastqc結(jié)束，如下；
運行MultiQC
直接指定MultiQC要分析的文件路徑即可，若數(shù)據(jù)在當前目錄下輸入multiqc .即可。
multiqc .
#使用--ignore忽略掉某些文件
#multiqc . --ignore *_R2*
multiqc結(jié)果解析，請學習：http://www.itdecent.cn/p/f83626fd1fa1
坑！坑！關(guān)于multiqc
我以為的又不是我以為的，2個小時出坑
一定要清楚multiqc是干啥的，我一直以為是和fastqc一樣的進行質(zhì)控，錯！我以為錯了，其實它是將fastqc質(zhì)檢后生成的.zip文件進行質(zhì)量信息結(jié)果合并呈現(xiàn)，所以我們應該在qc這個路徑下對.zip文件進行multiqc .
頓悟這個問題的心里路程：
當我第一次運行multiqc .最先出現(xiàn)這個報錯，下面繼續(xù)出錯顯示W(wǎng)ARNING
頓時蒙圈了，別人也沒這個情況呀，我就開始搜索，搜呀搜呀搜呀搜，第一個警告解決，幫助來自https://github.com/yaml/pyyaml/wiki/PyYAML-yaml.load(input)-Deprecation#footnotes
找到config.py文件，對其進行編輯，vi 路徑/config.py
在vi編輯器中查找yaml.load，一般都是yaml.load（f），這是我們要在其后面加上Loader=yaml.FullLoader，把所有的yaml.load中都加上，之后保存wq，退出vi
yaml.load (f, Loader=yaml.FullLoader)
接下來我們開始解決WARNING，用了1個小時解決的，我都想爆粗口，現(xiàn)在想想就一句話，multiqc處理的是fastqc質(zhì)控后生成的.zip，找對路徑找對文件執(zhí)行就完事了，要被自己蠢哭了，最后生成multiqc_data和.html，使用WinSCP打開.html
