這里是佳奧！經(jīng)歷了一個(gè)學(xué)期，做了幾遍轉(zhuǎn)錄組，感受頗深。
原理沒(méi)變，方法沒(méi)變，但是思維需要轉(zhuǎn)變，我們需要時(shí)刻保持清醒：
##這一步的代碼是做什么
##我當(dāng)前分析的是哪一個(gè)步驟
##這對(duì)于我的分析目的是否會(huì)有影響
做自己的項(xiàng)目不再是機(jī)械的重復(fù)步驟，而是要對(duì)自己的每一步負(fù)責(zé)，對(duì)課題組的經(jīng)費(fèi)支出負(fù)責(zé)。
那么，相信經(jīng)過(guò)了先前的學(xué)習(xí)，我們已經(jīng)掌握了轉(zhuǎn)錄組分析的全流程，這里，我會(huì)把使用的代碼全部整理出來(lái)，分成上下篇，希望對(duì)做非模式物種的同學(xué)們有所幫助。
那么我們開(kāi)始吧！
1 上傳實(shí)驗(yàn)數(shù)據(jù)
##首先在服務(wù)器的目錄下建立自己的目錄
/mnt/disk/lja/project/
mkdir tree
cd tree
##我們后續(xù)所有的分析都在tree目錄下
##上傳數(shù)據(jù)推薦使用MobaXterm左側(cè)文件欄的上傳（向上箭頭）功能
mkdir raw_fq
cd raw_fq
2 fastqc 質(zhì)量控制
##首先激活conda小環(huán)境
conda activate rnaseq
##批量質(zhì)控
ls *gz | xargs fastqc -t 10
##生成合成報(bào)告
multiqc ./
3 TrimGalore(依賴cutadapt) 過(guò)濾雙端測(cè)序結(jié)果
##安裝包下載，下載.gz至Linux下解壓即可
https://github.com/FelixKrueger/TrimGalore/releases
##解壓
tar -zxvf TrimGalore-0.6.7
##個(gè)人習(xí)慣，使用二進(jìn)制版，且每次都要添加到環(huán)境變量
export PATH="$PATH:/mnt/disk/lja/biosoft/trimgalory/TrimGalore-0.6.7"
##批處理過(guò)濾
##先生成config文件
ls *R1.fastq.gz >1
ls *R2.fastq.gz >2
paste 1 2 > config
##設(shè)置文件所在路徑
dir='/mnt/disk/lja/project/clean_fq'
##nohup掛起運(yùn)行，可top查看進(jìn)程情況
cat config | while read id
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
echo $dir $fq1 $fq2
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --paired -o $dir $fq1 $fq2 &
done
4 hisat2雙端比對(duì)(這一步前備份質(zhì)控后fq)
為什么這么說(shuō)呢？因?yàn)楸葘?duì)的時(shí)候很吃內(nèi)存和時(shí)間，中間要是有意外斷了的話，fq文件容易損壞，就無(wú)法再繼續(xù)分析了，就要回上一步重新過(guò)濾，為節(jié)省不必要的時(shí)間支出，盡量備份一遍。
讓我們繼續(xù)。
##到這一步，需要我們非模式種Tree的參考基因組文件，請(qǐng)找服務(wù)器管理員或者老師詢問(wèn)參考文件的位置
##建立索引，根據(jù)基因組大小，時(shí)間數(shù)小時(shí)不等
hisat2-build -p 8 Tree_V1.fasta genome
##原始文件名：Tree_1_val_1.fq.gz
##批量雙端比對(duì)
cd /clean_fq
ls *.gz | cut -d "_" -f 1 | sort -u | while read id ; do
ls -lh ${id}_1_val_1.fq.gz ${id}_2_val_2.fq.gz
nohup hisat2 -p 8 -x /mnt/disk/lja/refer/genome -1 ${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz -S ${id}.hisat2.sam
done
##批量轉(zhuǎn).bam
ls *.sam | while read id ; do (samtools sort -O bam -@ 5 -o $(basename ${id} ".sam").bam $id); done
##構(gòu)建.bam索引（如果.bam太大會(huì)無(wú)法建立）
ls *.bam | xargs -i samtools index {}
5 subread (featureCounts)
##上有分析最后一步，統(tǒng)計(jì)count數(shù)，得到rawcount矩陣
##這一步需要我們非模式種Tree的.gtf注釋文件
批量bam featureCounts:
gtf='/mnt/disk/lja/refer/Tree_V1.gtf'
nohup featureCounts -T 5 -p \
-a $gtf -o all.counts.txt \
/mnt/disk/lja/project/align/*.bam
至此，我們拿到了all.counts.txt，將它傳輸?shù)轿覀兊膫€(gè)人電腦，導(dǎo)入R開(kāi)始下游分析吧！
6 寫(xiě)在后面
上游分析相對(duì)比較簡(jiǎn)單，雖然涉及的軟件較多，但是思路是很明確的。
首先初步查看數(shù)據(jù)質(zhì)量，進(jìn)行過(guò)濾，再查看質(zhì)量，開(kāi)始與參考基因組比對(duì)，隨后count比對(duì)數(shù)值。
上述代碼適用于雙端測(cè)序結(jié)果，以及有參考基因組的情況。無(wú)參情況請(qǐng)翻閱先前轉(zhuǎn)錄組的博客。
那么，上游分析有什么值得注意的呢？
1 文件規(guī)范命名
一個(gè)規(guī)范的命名可以節(jié)省很多不必要找文件的時(shí)間。尤其在Linux系統(tǒng)下，找文件需要更加多的操作。
##一個(gè)新項(xiàng)目一個(gè)目錄
raw_fq：原始的fq文件
raw_qc：原始fq質(zhì)控報(bào)告
clean_fq：過(guò)濾后的fq文件
clean_qc：過(guò)濾后的質(zhì)控報(bào)告
align：比對(duì)（.sam .bam）
count：count矩陣
2 代碼報(bào)錯(cuò)一步步排查
首先，文件是否存在，目錄是否正確。
其次，conda環(huán)境是否正確激活，是否有該軟件，該軟件是否已被后臺(tái)占用 。
最后，--help查看軟件參數(shù)，是否是因?yàn)楦聦?dǎo)致的參數(shù)變動(dòng)。
3 固定流程
在沒(méi)有出現(xiàn)重大分析問(wèn)題的情況下，上游分析越快越好，因?yàn)槲覀兏枰褧r(shí)間花在下游分析上。
不僅是整理表達(dá)矩陣，還有在線注釋蛋白.fa數(shù)據(jù)，以及不斷篩選結(jié)果以對(duì)應(yīng)我們的實(shí)驗(yàn)變量。
那么，我們非模式種轉(zhuǎn)錄組下游分析篇再見(jiàn)！
