RNA-seq的上游分析主要基于Linux平臺(tái)，根據(jù)給定的fastq文件，經(jīng)過系列生信軟件處理最終得到表達(dá)矩陣的過程。原始數(shù)據(jù)文件一般可以從文獻(xiàn)中獲取（可能是sra格式，需要轉(zhuǎn)換成fastq，之后會(huì)專門整理下如何獲?。Ｎ抑熬毩?xí)的是老師直接給我的fastq文件(猴子的數(shù)據(jù))，還要補(bǔ)充的是雙末端測(cè)序同一樣本會(huì)有兩個(gè)文件的。就從這開始吧~

0、環(huán)境配置所需軟件

conda activate rna-seq
conda install fastqc -y  
conda install multiqc -y
conda install trimmomatic -y
conda install hisat2 -y
conda install trim-galore -y
conda install samtools -y
conda install -c bioconda subread

• 關(guān)于conda的說明，見之前的整理

1、質(zhì)控

目的：檢測(cè)fastq測(cè)序質(zhì)量如何，不好的話要進(jìn)行一定的處理才能用以分析
涉及軟件：fastqc、multiqc

ls *gz | xargs fastqc -o /home/monkey4/myfastqc -t 10

• 注意此時(shí)所在位置要在fastq文件路徑下
• xargs 命令實(shí)現(xiàn)了將ls的結(jié)果作為fastqc的參數(shù)輸入
• fastqc的 -o選項(xiàng)交代結(jié)果儲(chǔ)存路徑； -t 選項(xiàng)交代線程數(shù)
• 關(guān)于qc的結(jié)果，每個(gè)fastq文件都會(huì)有一個(gè)對(duì)應(yīng)的“檢查報(bào)告”，也可以將所有報(bào)告進(jìn)行合并。

cd /home/monkey4/myfastqc
multiqc ./
# multiqc會(huì)進(jìn)行自動(dòng)識(shí)別合并

• 關(guān)于報(bào)告結(jié)果解讀，網(wǎng)上有很多總結(jié)的介紹。之后自己可能也會(huì)整理下。目前的測(cè)序結(jié)果都蠻不錯(cuò)的，比如我做的那個(gè)。
• 如果的確測(cè)序結(jié)果質(zhì)量不行，可以進(jìn)行下面的處理。（我自己沒有實(shí)操過..僅記錄下）

#單獨(dú)處理
trim_glaore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o 輸出路徑 路徑/xxx_1.fastq.gz 路徑/xxx_2.fastq.gz

#批量處理
ls 路徑/*_1.fastq.gz >1
ls 路徑/*_2.fastq.gz >2
paste 1 2  #查看
paste 1 2 > config
cat >qc.sh
conda activate rna-seq
bin_trim_glore=trim_glaore
dir='/home/li/monkey4/test/clean'
cat $1 |while read id     #注意這里的$1 表示第一個(gè)參數(shù)，即bash qc.sh config 中的config
do
    arr=($id)
    fq1=${arr[0]}
    fq2=${arr[1]}
nohup $bin_trim_glore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 -paired -o $dir $fq1 $fq2 &
done
source deactivate
ctrl+c  #終止輸入編輯
bash qc.sh config

2、比對(duì)

• 目的：將fastq文件的所有read短序列比對(duì)到參考基因組上，即read本來屬于哪條基因上的。
• 涉及軟件：hisat2（還有同類型軟件bowtie2等）、samtools

2.1、構(gòu)建索引

構(gòu)建索引的原因是因?yàn)閰⒖蓟蚪M文件太太太大了，把一個(gè)幾百bp的read比對(duì)到基因組上，難度可想而知。而索引文件就類似圖書的目錄，按圖索驥就方便很多。

hisat2-build -p 10 基因組文件路徑/xxxx.fa 索引儲(chǔ)存路徑/xxxx
#索引文件命名不需要加后綴，取個(gè)名即可。軟件會(huì)自己加的
hisat2-build -p 4 /home/li/genome/monkey.fa  /home/li/monkey4/hisat2-index/Macaca
#這個(gè)是我的實(shí)操代碼

2.2、進(jìn)行比對(duì)

#單個(gè)比對(duì)
hisat2 -t -p 4 -x 索引文件路徑/索引文件名 -1 XXX1.fq.gz -2 XXX2.fq.gz -S 儲(chǔ)存路徑/XXX.sam
# 針對(duì)雙末端測(cè)序，可以同時(shí)比對(duì)，進(jìn)行結(jié)果合并
hisat2 -t -p 4 -x /home/li/monkey4/rnaseq_test/index/hisat2_index/Macaca -1 M5Acortex_1.fq.gz -2 M5Acortex_2.fq.gz -S /home/li/monkey4/rnaseq_test/M5A.sam
#實(shí)操代碼

#批量比對(duì)
ls *gz |cut -d"_" -f 1 | sort -u |while read id ; do hisat2 -p 10 -x 索引文件路徑/索引文件名 -1 ${id}_1.fq.gz -2 ${id}_2.fq.gz  -S 輸出路徑/${id}.hisat.sam ; done
ls *gz |cut -d"_" -f 1 | sort -u |while read id ; do hisat2 -p 10 -x /home/li/monkey4/rnaseq_test/index/hisat2_index/Macaca -1 ${id}_1.fq.gz -2 ${id}_2.fq.gz -S /home/li/monkey4/rnaseq_test/sam/${id}.sam 1>$id.download.log 2>&1 ; done

分析得到的sam文件記錄一個(gè)樣本的所有read比對(duì)信息。

2.3、sam轉(zhuǎn)bam

samtools view -S XXX.sam -@ 10 -b > XXX.bam
ls *.sam | cut -d"." -f 1 | sort -u | while read id ;do samtools view -S ${id}.hisat.sam -@ 10 -b > ${id}.bam ; done

而bam格式中的b是binary的意思，是sam格式的二進(jìn)制表示方式，主要就是為了節(jié)省內(nèi)存，提高后面分析速度。

3、計(jì)數(shù)

• 目的：統(tǒng)計(jì)每個(gè)基因比對(duì)上了多少個(gè)read
• 涉及軟件：subread

#單個(gè)計(jì)數(shù)
featureCounts -T 10 -p -t exon -g gene_id -a 基因組注釋文件路徑/XXX.gtf -o 結(jié)果文件名 bam文件路徑/XXX.bam
featureCounts -T 10 -p -t exon -g gene_id -a /home/li/genome/Macaca_fascicularis.Macaca_fascicularis_5.0.97.gtf -o counts.M2Ccortex.txt /home/li/monkey4/test/aligned/M2Ccortex.bam

#批量處理
ls *tex.bam |cut -d"." -f 1 |sort -u |while read id ; do featureCounts -T 10 -p -t exon -g gene_id -a /home/li/genome/Macaca_fascicularis.Macaca_fascicularis_5.0.97.gtf -o counts.${id}.txt /home/li/monkey4/test/aligned/${id}.bam; done

#將批量計(jì)數(shù)結(jié)果轉(zhuǎn)成CSV格式
ls *.txt |cut -d"." -f 1,2 |sort -u | while read id ; do cut -f 1,7-12 ${id}.txt | grep -v "#" > countmatrix.${id}.csv; done

#合并所有樣本的CSV文件
paste *.csv | awk '{printf $1 "\t" ; for (i=2; i<=NF; i=i+2) printf $i "\t" ; print $i}' > merge.count.csv

以上就是一個(gè)極簡(jiǎn)版的RNA-seq的上游分析過程，最終獲得了基因表達(dá)矩陣，用以后期的數(shù)據(jù)分析。文中提到一些數(shù)據(jù)文件及格式還是值得后期慢慢梳理的。暫時(shí)先這樣吧。上述代碼主要參考自生信技能樹的學(xué)習(xí)，站在巨人肩膀?qū)W生信，嘻嘻~