寫在開頭：在這一個合集中，我將詳細(xì)介紹bulk RNA-seq實(shí)際數(shù)據(jù)分析的流程闻伶，有哪些注意的地方以及一些小技巧

首先是獲取測序數(shù)據(jù)：如果是自己的測序數(shù)據(jù)，直接從公司給的賬號中下載各個樣本的RawData（fastq格式）即可毅桃。
如果想要分析公共數(shù)據(jù)珍手，在GEO數(shù)據(jù)庫輸入文章中給出的登錄號，如GSE184771双抽，這些數(shù)據(jù)通常是經(jīng)過質(zhì)控百框、數(shù)據(jù)歸一化、差異表達(dá)分析等處理牍汹。提交者通常會詳細(xì)描述其數(shù)據(jù)的處理步驟铐维。Series Matrix File可能會提供關(guān)于實(shí)驗(yàn)設(shè)計(jì)、數(shù)據(jù)處理的詳細(xì)信息慎菲；supplementary file可能會提供經(jīng)過處理的表達(dá)矩陣嫁蛇，也就是我們需要下載到R中進(jìn)行分析的文件÷陡茫或者進(jìn)入下面的SRA Run selector睬棚，SRA數(shù)據(jù)庫一般儲存的是原始的測序文件，下載SRR_Acc_List1.txt解幼，以便在Linux中下載SRA數(shù)據(jù)抑党，可以使用它們進(jìn)行完整的上游和下游分析。

GEO184771

Accession list

#SRAToolkit 是NCBI提供的用于處理來自SRA 數(shù)據(jù)庫測序數(shù)據(jù)的一個工具包撵摆。前面的常用Linux命令合集中已經(jīng)介紹如何下載了
#使用其中的prefetch命令批量下載數(shù)據(jù)
#運(yùn)行nohup  &  用于放入后臺下載底靠，避免關(guān)閉終端導(dǎo)致下載中斷
mkdir -p GSE184771/fastq
cd /home/GSE184771/fastq
#讀取文件 SRR_Acc_List1.txt 的內(nèi)容，$()將讀取到的內(nèi)容作為參數(shù)傳遞到prefetch命令中特铝，最后-O .指定下載文件的輸出目錄為當(dāng)前目錄
nohup prefetch -O . $(<SRR_Acc_List1.txt) &

#下載完 sra 文件暑中，我們往往需要先將 SRA文件轉(zhuǎn)化為 fastq 文件才能用于后續(xù)分析。使用前面安裝的 SRAToolkit 中的fastq-dump工具
vi RNA-seq.sh
# 讀取包含本地 SRA 文件路徑的文件SRR_Acc_List2.txt苟呐，并逐行處理
#其中fastq-dump參數(shù)：--split-files將配對reads（paired-end reads）分別輸出到不同的文件中痒芝，對于一方有而一方?jīng)]有的reads直接丟棄俐筋。--split-3將配對reads輸出到兩個文件中牵素，并將unpaired reads輸出到第三個文件中。--split-spot將雙端測序分為兩份,但是都放在同一個文件中澄者。
mkdir /home/GSE184771/1.QC
cd /home/GSE184771/fastq
cat SRR_Acc_List2.txt | while read i
do
    # fastq-dump 命令處理SRA 文件笆呆，并將其轉(zhuǎn)換為壓縮的 FASTQ 文件，輸出到當(dāng)前目錄粱挡，fastq-dump 命令處理本地的 SRA 文件赠幕，并將其轉(zhuǎn)換為壓縮的 FASTQ 文件，輸出到當(dāng)前目錄
    fastq-dump --gzip --split-files $i -O /home/GSE184771/1.QC && echo "** ${i} to fastq done **"
done
bash RNA-seq.sh
fastq-dump --split-files ./SRR16060501/SRR16060501.sra -O /home/GSE184771/1.QC #或者一個一個運(yùn)行不報(bào)錯

現(xiàn)在均得到了樣本的原始測序fastq文件询筏，接下來進(jìn)行質(zhì)控榕堰。

#如果樣本數(shù)不多，其實(shí)可以一個一個運(yùn)行
#列出所有（_R1.fq.gz）的路徑，并將這些路徑保存到文件1中逆屡。>符號被用作輸出重定向操作符圾旨，作用是將命令的輸出結(jié)果寫入到指定的文件中。
ls /home/RNA-seq/fastq/*_R1.fq.gz > 1 
ls /home/RNA-seq/fastq/*_R2.fq.gz > 2
#使用cut命令根據(jù)/分隔符提取第5個字段（第一個字段為空魏蔗，完整文件路徑在第5個位置）砍的，再次使用cut根據(jù)_分隔符提取第1個字段（樣本名），并將結(jié)果保存到文件0中莺治。
ls /home/RNA-seq/fastq/*_R2.fq.gz |cut -d"/" -f 5|cut -d"_" -f 1  > 0
#paste命令合并文件0廓鞠、1和2的內(nèi)容（分別是樣本名、_R1.fq.gz和_R2.fq.gz的路徑）谣旁，生成文件config
paste 0 1 2  > config
#類似于SRR16060501 /home/RNA-seq/fastq/SRR16060501_R1.fq.gz /home/RNA-seq/fastq/SRR16060501_R2.fq.gz
vi RNA-seq.sh
for i in $(cat config)
do
    echo $i
    arr=($i)
    fq1=${arr[1]}
    fq2=${arr[2]}
    sample=${arr[0]}
      fastqc fq1 fq2 -o /home/GSE184771/1.QC
done
cd /home/GSE184771/1.QC
#整合質(zhì)控結(jié)果
multiqc ./

得到了質(zhì)控結(jié)果床佳，我以一個樣本的R1為例。QC結(jié)果包括以下幾部分榄审。

QC結(jié)果

其中比較重要的包括：per base sequence quality每個堿基的質(zhì)量夕土。藍(lán)色的線表示每個堿基位置上的平均質(zhì)量分?jǐn)?shù)（mean quality score），一般還會有黃色的框（箱線圖）瘟判，表示的是每個堿基位置的質(zhì)量分?jǐn)?shù)分布的中位數(shù)和四分位范圍（interquartile range, IQR）怨绣，一般要求藍(lán)線要在綠色區(qū)域內(nèi)，箱線圖在黃色區(qū)域以上拷获。堿基質(zhì)量分?jǐn)?shù)與堿基出錯的概率關(guān)系為Q=?10log10(P)篮撑，某個位置的堿基質(zhì)量為30就是它出錯的概率為0.1%

per base sequence quality

Per base sequence content

Sequence duplication levels峰一般在最左側(cè)赢笨。圖中顯示有50%左右的reads有大于10個重復(fù)。測序數(shù)據(jù)中存在大量的重復(fù)序列驮吱，有可是PCR擴(kuò)增偏好茧妒；測序深度不足，可能會導(dǎo)致樣本中的某些序列被過度重復(fù)(當(dāng)測序深度不足時左冬，測序過程中的隨機(jī)抽樣誤差會更顯著桐筏。測序儀在讀取DNA片段時是隨機(jī)的，如果測序深度低拇砰，每個位置被讀取的次數(shù)有限梅忌，這會導(dǎo)致某些片段被重復(fù)讀取，而其他片段可能完全沒有被讀取除破。)牧氮；生物學(xué)因素，某些生物樣本中天然存在高重復(fù)率的序列瑰枫，特別是在低復(fù)雜度區(qū)域或高表達(dá)的基因區(qū)域踱葛；技術(shù)問題如接頭污染或樣本準(zhǔn)備中的錯誤，也可能導(dǎo)致重復(fù)序列的增加。

Sequence duplication levels

Overrepresented sequences顯示這條序列大量出現(xiàn)尸诽，占全部reads的13%左右圾笨，這很可能就是接頭，或者PCR偏好性擴(kuò)增造成的逊谋，一般后續(xù)需要過濾掉

Overrepresented sequences

Adapter content擂达，如果有接頭殘留，后續(xù)一定要過濾掉

adapter content

接下來就需要針對質(zhì)控情況進(jìn)行選擇性過濾了胶滋，使用Trim-galore對數(shù)據(jù)進(jìn)行過濾和質(zhì)控板鬓，它將fastqc與cutadapt封裝在一起。

#這個可以針對不同文件質(zhì)控情況究恤，選擇不同的過濾參數(shù)進(jìn)行過濾
mkdir /home/GSE184771/2.trim
cd /home/GSE184771/fastq
#使用--paired 參數(shù)時俭令，--adapter 參數(shù)會應(yīng)用于R1，--adapter2 參數(shù)會應(yīng)用于R2部宿，從而去除自定義的序列抄腔，比如overrepresented sequences。--illumina 參數(shù)會同時處理R1和R2文件中的Illumina接頭序列理张。--quality設(shè)置質(zhì)量分?jǐn)?shù)閾值赫蛇，低于該閾值的堿基將被修剪。--stringency 4設(shè)置修剪的嚴(yán)格度雾叭，至少有 4 個匹配adapter的堿基序列才能進(jìn)行修剪悟耘。-j 8使用的線程數(shù)為 8。-o 輸出目錄织狐。--fastqc會對修剪后的文件運(yùn)行 FastQC 分析
nohup trim_galore --paired --illumina --adapter2 AAGCAGTGGTATCAACGCAGAGTACTTTTTTTTTTTTTTTTTTTTTTTTT --fastqc C1_R1.fq.gz C2_R2.fq.gz --quality 25 --stringency 4 -j 8 --fastqc -o /home/GSE184771/2.trim &
#其他根據(jù)情況進(jìn)行選擇性過濾
#生成的_R1_val_1.fq.gz和_R2_val_2.fq.gz: 表示經(jīng)過修剪并驗(yàn)證配對關(guān)系后的R1和R2文件,是用于后續(xù)比對和分析的文件

Bowtie2 和 HISAT2 是用于高通量測序數(shù)據(jù)的序列比對工具暂幼。Bowtie2對短讀長數(shù)據(jù)（50-100bp）比對非常高效，支持局部（local）和全局（global）比對模式移迫，可以用于ChIP-seq旺嬉、DNA-seq以及小 RNA-seq。HISAT2 專門用于處理 RNA-seq 數(shù)據(jù)厨埋；適合長讀長數(shù)據(jù)邪媳，優(yōu)化了比對長讀長（例如 100-300bp）的能力，同時也可以高效處理基因組中的剪接事件揽咕。所以我們這里使用Hisat2進(jìn)行比對

#需要準(zhǔn)備參考物種的參考基因組和注釋文件悲酷，有很多網(wǎng)站上的數(shù)據(jù)可供下載套菜，如UCSC亲善、GENCODE、Ensembl逗柴、NCBI蛹头、iGenomes等。對于大多數(shù)基因組比對和下游分析，primary assembly的 FASTA和注釋 文件通常已經(jīng)足夠
#首先構(gòu)建索引渣蜗，輸入?yún)⒖蓟蚪M序列屠尊，構(gòu)建的索引的基礎(chǔ)名稱為index_m39.genome。時間長建議后臺運(yùn)行
mkdir /home/GSE184771/reference
mkdir /home/GSE184771/3.align
cd /home/GSE184771/reference
#在構(gòu)建 HISAT2 索引時耕拷，可以使用 --ss 和 --exon 選項(xiàng)來包含剪接位點(diǎn)和外顯子信息
#使用hisat2_extract_splice_sites.py 和 hisat2_extract_exons.py 腳本從 GTF 文件中提取剪接位點(diǎn)和外顯子信息讼昆。
hisat2_extract_splice_sites.py gencode.GRCm39_primary.annotation.gtf > m39.ss
hisat2_extract_exons.py gencode.GRCm39_primary.annotation.gtf > m39.exon
nohup hisat2-build -p 4 --ss m39.ss --exon m39.exon GRCm39.primary_genecode.genome.fa index_m39.genome &
#-x為構(gòu)建的索引。-1和-2為過濾后的文件骚烧，-S為輸出路徑浸赫，--dta啟用下游轉(zhuǎn)錄組分析，適用于轉(zhuǎn)錄組數(shù)據(jù)赃绊。--new-summary生成詳細(xì)的比對統(tǒng)計(jì)信息
nohup hisat2 -x /home/GSE184771/reference/index_m39.genome/index_m39.genome -1 C1_R1_val_1.fq.gz -2 C1_R2_val_2.fq.gz -S /home/GSE184771/3.align/C1_aligned.sam -p 8 --dta --new-summary &
#將比對結(jié)果壓縮成bam文件并根據(jù)染色體排序既峡，以節(jié)省存儲空間、提高處理效率碧查。排序后的 BAM 文件可以被索引运敢，索引后的 BAM 文件可以快速隨機(jī)訪問特定區(qū)域的數(shù)據(jù)，這對于下游的分析和可視化非常重要忠售，很多下游分析工具要求輸入的 BAM 文件是按坐標(biāo)排序的（如stringtie)
cd /home/GSE184771/3.align
nohup samtools sort -@ 8 C1_aligned.sam -o C1_sorted.bam &
建立BAM索引文件
nohup samtools index ./C1_sorted.bam &

reference fasta

Stringtie定量传惠，gtf文件不能是壓縮文件

mkdir /home/GSE184771/4.stringtie
cd /home/GSE184771
#-G指定GTF文件，-e 參數(shù)表示只對提供的注釋文件中的轉(zhuǎn)錄本進(jìn)行表達(dá)量估算稻扬，而不進(jìn)行新的轉(zhuǎn)錄本組裝涉枫。C1_genes.list 文件是一個基因級別的表達(dá)量表格，列出了每個基因的基本信息和表達(dá)量信息腐螟。C1_count.gtf 文件是一個轉(zhuǎn)錄本級別的 GTF 文件愿汰，包含了轉(zhuǎn)錄本和外顯子的信息及其表達(dá)量。
nohup stringtie ./3.align/C1_sorted.bam -G ./reference/gencode.GRCm39.annotation.gtf -e -p 8 -o ./4.stringtie/C1_count.gtf -A ./stringtie/C1_genes.list &
#由于后續(xù)R下游分析使用到的基因表達(dá)矩陣都是原始的count矩陣乐纸，不能是經(jīng)過FPKM或TPM歸一化衬廷，以及其他轉(zhuǎn)換后的數(shù)據(jù)。所以需要從中得到原始基因表達(dá)矩陣
#匯總stringtie表達(dá)量結(jié)果汽绢，例如
vi sample_list.txt
C1  /home/zhangyina/zbw_RNA-seq/stringtie/C1_count.gtf
C2  /home/zhangyina/zbw_RNA-seq/stringtie/C2_count.gtf
C3  /home/zhangyina/zbw_RNA-seq/stringtie/C3_count.gtf
E1  /home/zhangyina/zbw_RNA-seq/stringtie/E1_count.gtf
E2  /home/zhangyina/zbw_RNA-seq/stringtie/E2_count.gtf
E3  /home/zhangyina/zbw_RNA-seq/stringtie/E3_count.gtf
#將如上內(nèi)容寫入sample_lsit.txt,該文件為 \t（Tab） 分隔的兩列吗跋，第一列為樣本名稱，第二列為定量的 GTF 文件的路徑宁昭，注意最后不要有空行
cat sample_lsit.txt
#從https://github.com/gpertea/stringtie跌宛，下載prepDE.py3腳本
python prepDE.py3 -i sample_list.txt
#執(zhí)行完上述步驟后得到的gene_count_matrix.csv即為基因表達(dá)矩陣，可以導(dǎo)入到R語言用edgeR / DESeq2包進(jìn)行差異表達(dá)基因分析积仗〗校基因表達(dá)矩陣包含了每個基因在不同樣本中的原始計(jì)數(shù)（即每個基因在每個樣本中檢測到的讀數(shù)數(shù)量）。行表示基因寂曹，列表示樣本哎迄。
#生成的transcript_count_matrix.csv 文件包含了每個轉(zhuǎn)錄本在不同樣本中的原始計(jì)數(shù)（即每個轉(zhuǎn)錄本在每個樣本中檢測到的讀數(shù)數(shù)量）回右。行表示轉(zhuǎn)錄本，列表示樣本漱挚。