今天開始正式進(jìn)入RNA-seq分析啦~
軟件安裝
由于昨天就已經(jīng)創(chuàng)建好了新的conda的環(huán)境摩渺，所以今天只要把用到的軟件都裝一下就好了贡这。
把要安裝的軟件分一下類型：
數(shù)據(jù)格式轉(zhuǎn)換
-
sra-tools: 這個(gè)軟件用于把NCBI里下載的.sra格式的數(shù)據(jù)給轉(zhuǎn)換成fastq赊瞬。用到的工具是里面的fastq-dump
數(shù)據(jù)質(zhì)控
數(shù)據(jù)質(zhì)量的好壞是數(shù)據(jù)分析中的關(guān)鍵，所謂“garbage in, garbage out”贮泞，如果一開始的數(shù)據(jù)是有問題的移迫，后續(xù)再怎么分析都是徒勞。
那數(shù)據(jù)會(huì)存在什么樣的問題呢败潦？
主要是二代測(cè)序增加進(jìn)去adapters序列本冲，和測(cè)序中的低質(zhì)量數(shù)據(jù)。它們會(huì)影響后續(xù)的比對(duì)劫扒、定量準(zhǔn)確性檬洞，更會(huì)影響下游分析，因此要先生成質(zhì)量報(bào)告了解一下數(shù)據(jù)的質(zhì)量沟饥，如果質(zhì)量不好的話要把不符合要求的序列給過濾掉添怔。
所以這一塊分兩步，1. 生成質(zhì)量報(bào)告贤旷；2. 過濾不符合要求的序列广料。
-
fastqc&multiqc: 用于生成質(zhì)量報(bào)告。fastqc能對(duì)每個(gè)樣品生成單獨(dú)的質(zhì)量報(bào)告幼驶，而multiqc則是可以把多個(gè)樣本的fastqc結(jié)果整合到一個(gè)文件中艾杏，這樣比較方便全面的了解和對(duì)比多個(gè)數(shù)據(jù)之間的質(zhì)量情況。 -
fastp: 由陳實(shí)富大佬寫的新手友好型質(zhì)控軟件盅藻，只要設(shè)置輸入文件和輸出文件就好了糜颠，推薦給各位胖友汹族。之前寫過一篇關(guān)于這個(gè)軟件的簡書文章：用fastp對(duì)轉(zhuǎn)錄組數(shù)據(jù)做QC -
trim_galore&trimmomatic另外兩個(gè)用于做質(zhì)控的軟件。為了防止質(zhì)控軟件本身的算法原因?qū)е路治鼋Y(jié)果的誤差其兴，條件允許的情況下建議用不同的質(zhì)控軟件來處理數(shù)據(jù)進(jìn)行對(duì)比顶瞒，取結(jié)果好的進(jìn)行下游分析。
比對(duì)（mapping）
有參轉(zhuǎn)錄組和無參轉(zhuǎn)錄組分析的本質(zhì)區(qū)別就在于有沒有這個(gè)“參”元旬。參是什么呢榴徐？是reference，中文世界一般叫“參考序列”匀归。而我研究生期間做的植物基因組坑资，就是在做這個(gè)“參”。
比對(duì)的軟件有很多種穆端，有基于基因組比對(duì)（star袱贮、hisat2）、基于轉(zhuǎn)錄組比對(duì)(bowtie体啰、bwa)攒巍、不基于比對(duì)(salmon)
我目前用bwa比較多，別的只是知道荒勇，沒怎么用過柒莉。
定量
-
subread用的是其下的featureCounts這個(gè)小程序。這一塊我不太了解沽翔，這也是我參加這次轉(zhuǎn)錄組學(xué)習(xí)小組的原因~
用conda安裝軟件
# 激活專門用于RNA-seq分析的conda環(huán)境
conda activate rnaseq
# 安裝軟件
conda install fastp fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon
因?yàn)槲夷壳叭嗽赟ingapore兢孝，所以用conda下載軟件還挺快的。有的小朋友習(xí)慣在conda install的時(shí)候加-y參數(shù)仅偎，這樣就只要等著它裝好就好了跨蟹，不需要再手動(dòng)去確定是否安裝這些軟件。但是也有個(gè)問題橘沥，它在整個(gè)安裝過程中不會(huì)給任何的提示和輸出窗轩，就只能看到solving environment 的狀態(tài)欄在“愛的魔力轉(zhuǎn)圈圈”，當(dāng)安裝很多軟件或者網(wǎng)絡(luò)不好的時(shí)候會(huì)特別沒有安全感威恼。
數(shù)據(jù)格式轉(zhuǎn)換
我想想我一般是怎么寫這個(gè)的：
#!/usr/bin/env bash
list="SRR1039509.sra
SRR1039512.sra
SRR1039513.sra
SRR1039516.sra
SRR1039517.sra
SRR1039520.sra
SRR1039521.sra"
for i in $list
do
fastq-dump --gzip --split-3 ${i} -O ../01raw
done
我的list是通過下面的命令生成的:
ls *.sra | xargs -l
而豆豆提供的腳本是這個(gè)畫風(fēng)：
raw=~/rnaseq/raw # 這個(gè)路徑根據(jù)個(gè)人需求進(jìn)行修改即可
cat $raw/srr.ids | while read i ;do
echo $i
# time fastq-dump --gzip --split-3 -A $i $raw/${i}.sra -O $raw;
done
確實(shí)看起來比我高明得多呢品姓。學(xué)習(xí)了學(xué)習(xí)了寝并。用time這個(gè)命令計(jì)算運(yùn)行時(shí)間是一個(gè)很好的習(xí)慣箫措。
先用井號(hào)注釋掉了真正要運(yùn)行的命令，這在寫腳本的過程中是個(gè)好習(xí)慣衬潦，因?yàn)橛袝r(shí)自己會(huì)搞錯(cuò)變量(比如這里的i) 斤蔓，于是用echo $i可以先看看是不是真正要用的，先運(yùn)行一下看看結(jié)果
把腳本寫進(jìn)一個(gè)叫做fastqdump.sh的文件镀岛，然后nohup運(yùn)行弦牡。
(不要隨便起個(gè)名字噢友驮，不然過段時(shí)間就忘了這個(gè)當(dāng)時(shí)寫了是做什么用的了呢，后患無窮驾锰，整理起來特別頭疼卸留。)
我用nohup主要是它可以保留軟件運(yùn)行過程中的日志和報(bào)錯(cuò)，當(dāng)然不用nohup也可以做到這點(diǎn)：
bash fastqdump.sh 1>fastqdump.log 2>&1
# 或者
bash fastqdump.sh 1>fastqdump.log 2>fastqdump.err
1>fastqdump.log 表示將結(jié)果的正確日志文件輸出到fastqdump.log中椭豫，2>&1 這個(gè)2表示錯(cuò)誤日志耻瑟，將它也合并到1的正確日志中。
當(dāng)然赏酥，如果不想要錯(cuò)誤信息喳整，可以把它丟到linux"黑洞"中去: 2>/dev/null
下載參考序列
根據(jù)文章介紹，所用的reference是hg19的版本裸扶。因此就下載這個(gè)版本的genome和gtf數(shù)據(jù)備用框都。
wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
太晚了。明天繼續(xù)寫呵晨。
Hans Chen
2019年6月9日于 Westwood Hostel, Jurong West Ave. 5, Singapore.
