本文參考文獻(xiàn)https://www.nature.com/articles/nprot.2016.095#materials進(jìn)行bulk RNA-seq數(shù)據(jù)分析

首先進(jìn)行原始數(shù)據(jù)的下載，當(dāng)然因為文章中已經(jīng)有了鏈接進(jìn)行數(shù)據(jù)的下載，因此直接使用wget命令就能完成所需數(shù)據(jù)的下載，但是也要學(xué)會基本的測序數(shù)據(jù)、參考基因組、注釋文件、index索引文件的下載，當(dāng)然也可以自己建立index文件。

raw data下載（找了2個小鼠的測序結(jié)果）：

fastq-dump --split-3 --gzip SRR3589959

fastq-dump --split-3 --gzip SRR3589960

加上--split-3之后, 會把原來雙端拆分成兩個文件,但是原來單端并不會保存成兩個文件. 還有你用--gzip就能輸出gz格式, 能夠節(jié)省空間的同時也不會給后續(xù)比對軟件造成壓力,比對軟件都支持，就是時間要多一點。

標(biāo)準(zhǔn)基因組的數(shù)據(jù)下載：

wgetftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.fna.gz

注釋文件的下載：

Wgetftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/635/GCF_000001635.26_GRCm38.p6/GCF_000001635.26_GRCm38.p6_genomic.gff.gz

Index的下載，自己的電腦沒法建立index有現(xiàn)成的index可以去https://ccb.jhu.edu/software/hisat2/index.shtml下載：

關(guān)于HISAT2官網(wǎng)上的index：https://www.biostars.org/p/290721/

genome: HFM index for reference

genome_snp: HGFM index for reference plus SNPs

genome_tran: HGFM index for reference plus transcripts

genome_snp_tran: HGFM index for reference plus SNPs and transcripts

當(dāng)然也可以自己建立索引文件：

Hisat2-build reference genome

Index文件不會用不就狗帶了嗎，注意index文件的格式，其實就是只有電腦能看你不能看的二進(jìn)制，所以千萬不要嘗試打開，否則會是一堆亂碼，感覺自己弄錯了emmm，打開index文件是一些.ht2格式的文件，參考基因組有幾條染色體就有幾個這樣的文件，但是運行HISAT2的時候千萬不要在文件名后面加上后綴！要不然會報錯的~

分析步驟：

[if !supportLists]1、[endif]首先是對測序數(shù)據(jù)質(zhì)量的檢測——fastqc

fastqc-o ~/rnaseq/data/firstqc -t 6 ERR188044_chrX_1.fastq.gz && fastqc -o~/rnaseq/data/firstqc -t 6 ERR188044_chrX_2.fastq.gz

2、去除接頭和低質(zhì)量reads

trim_galore-output_dir ~/rnaseq/data/firstqc/clean --paired --length 58 --quality 25ERR188044_chrX_1.fastq.gz ERR188044_chrX_2.fastq.gz

[if !supportLists]2、[endif]進(jìn)行第二次質(zhì)檢：

fastqc-o ~/rnaseq/data/secondqc -t 6 ERR188044_chrX_1_val_1.fq.gzERR188044_chrX_2_val_2.fq.gz

結(jié)果和原來差不多，因此便用了原來的reads進(jìn)行比對。

[if !supportLists]3、[endif]將reads比對到參考基因組

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188044_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188044_chrX_2.fastq.gz -S ERR188044_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188104_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188104_chrX_2.fastq.gz -S ERR188104_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188234_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188234_chrX_2.fastq.gz -S ERR188234_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188245_chrX_1.fastq.gz -2 ~/rnaseq/chrX_data/samples/ERR188245_chrX_2.fastq.gz-S ERR188245_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188257_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188257_chrX_2.fastq.gz -S ERR188257_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188273_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188273_chrX_2.fastq.gz -S ERR188273_chrX.sam

hisat2 -p 6 --dta -t -x ~/rnaseq/chrX_data/indexes/chrX_tran-1 ~/rnaseq/chrX_data/samples/ERR188337_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188337_chrX_2.fastq.gz -S ERR188337_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188383_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188383_chrX_2.fastq.gz -S ERR188383_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188401_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188401_chrX_2.fastq.gz -S ERR188401_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1~/rnaseq/chrX_data/samples/ERR188428_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188428_chrX_2.fastq.gz -S ERR188428_chrX.sam

hisat2 -p 6 --dta -t -x ~/rnaseq/chrX_data/indexes/chrX_tran-1 ~/rnaseq/chrX_data/samples/ERR188454_chrX_1.fastq.gz -2~/rnaseq/chrX_data/samples/ERR188454_chrX_2.fastq.gz -S ERR188454_chrX.sam

hisat2 -p 6 --dta -t -x~/rnaseq/chrX_data/indexes/chrX_tran -1 ~/rnaseq/chrX_data/samples/ERR204916_chrX_1.fastq.gz-2 ~/rnaseq/chrX_data/samples/ERR204916_chrX_2.fastq.gz -S ERR204916_chrX.sam

以上參數(shù)都可以通過hisat2 –help查到。

4、將sam文件轉(zhuǎn)換成bam文件并進(jìn)行排序，至于為什么要進(jìn)行排序是因為打開其中一個sam文件就可以看到人家已經(jīng)指明數(shù)據(jù)是unsorted的了。

samtoolssort -@ 8 -o ~/bam/ERR188044_chrX.bam ERR188044_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188104_chrX.bam ERR188104_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188234_chrX.bam ERR188234_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188245_chrX.bam ERR188245_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188257_chrX.bam ERR188257_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188273_chrX.bam ERR188273_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188337_chrX.bam ERR188337_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188383_chrX.bam ERR188383_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188401_chrX.bam ERR188401_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188428_chrX.bam ERR188428_chrX.sam && samtoolssort -@ 8 -o ~/bam/ERR188454_chrX.bam ERR188454_chrX.sam

[if !supportLists]4、[endif]進(jìn)行轉(zhuǎn)錄本的組裝。

1092? stringtie -p 8 -G chrX.gtf -oERR188044_chrX.gtf -l ERR188044 ERR188044_chrX.bam

?1093? ls

?1094?stringtie -p 8 -G chrX.gtf -o ERR188104_chrX.gtf -l ERR188104ERR188104_chrX.bam

?1095? ls

?1096?stringtie -p 8 -G chrX.gtf -o ERR188234_chrX.gtf -l ERR188234ERR188234_chrX.bam

?1097? ls

?1098?stringtie -p 8 -G chrX.gtf -o ERR188245_chrX.gtf -l ERR188245ERR188245_chrX.bam

?1099?stringtie -p 8 -G chrX.gtf -o ERR188257_chrX.gtf -l ERR188257ERR188257_chrX.bam

?1100? ls

?1101?stringtie -p 8 -G chrX.gtf -o ERR188273_chrX.gtf -l ERR188273ERR188273_chrX.bam

?1102?stringtie -p 8 -G chrX.gtf -o ERR188337_chrX.gtf -l ERR188337ERR188337_chrX.bam

?1103? ls

?1104?stringtie -p 8 -G chrX.gtf -oERR188383_chrX.gtf -l ERR188383 ERR188383_chrX.bam

?1105?stringtie -p 8 -G chrX.gtf -o ERR188401_chrX.gtf -l ERR188401ERR188401_chrX.bam

?1106?stringtie -p 8 -G chrX.gtf -o ERR188428_chrX.gtf -l ERR188428ERR188428_chrX.bam

?1107?stringtie -p 8 -G chrX.gtf -o ERR188454_chrX.gtf -l ERR188454ERR188454_chrX.bam

[if !supportLists]5、[endif]將組裝后的轉(zhuǎn)錄本merge一下：

stringtie --merge -p8 -G chrX.gtf -o stringtie_merged.gtf ~/rnaseq/chrX_data/mergelist.txt

其中mergelist.txt就是上一步得到的文件名字組成的純文本格式的文件，當(dāng)然也可以將上述文件名直接打上去，以空格隔開。下圖為mergelist.txt文件打開的樣子：

[if !supportLists]6、[endif]將merge后的注釋文件與標(biāo)準(zhǔn)基因組的注釋文件進(jìn)行對比（這一步可以不做）：

gffcompare -rchrX.gtf -G -o merged stringtie_merged.gtf

[if !supportLists]7、[endif]再以merge后的注釋文件作為參考進(jìn)行轉(zhuǎn)錄本的組裝并輸出為ballgown可以識別的格式（因為在一個實驗中，我們往往有很多樣本，而不同樣本的轉(zhuǎn)錄本往往是不同的，而我們需要讓他們在近似相同的情況下進(jìn)行組裝，因此再用merge后的注釋文件作為參考進(jìn)行組裝）：

1117? stringtie–e –B -p 8 -G stringtie_merged.gtf -o ~/ballgown/ERR188044/ERR188044_chrX.gtfERR188044_chrX.bam

?1118? stringtie -e –B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam

?1119? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188044/ERR188044_chrX.gtf ERR188044_chrX.bam

?1120? ls

?1121? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188104/ERR188104_chrX.gtf ERR188104_chrX.bam

?1122? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188234/ERR188234_chrX.gtf ERR188234_chrX.bam

?1123? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188245/ERR188245_chrX.gtf ERR188245_chrX.bam

?1124? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188257/ERR188257_chrX.gtf ERR188257_chrX.bam

?1125? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188273/ERR188273_chrX.gtf ERR188273_chrX.bam

?1126? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188337/ERR188337_chrX.gtf ERR188337_chrX.bam

?1127? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188383/ERR188383_chrX.gtf ERR188383_chrX.bam

?1128? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188401/ERR188401_chrX.gtf ERR188401_chrX.bam

?1129? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188428/ERR188428_chrX.gtf ERR188428_chrX.bam

?1130? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR188454/ERR188454_chrX.gtf ERR188454_chrX.bam

?1131? stringtie -e -B -p 8 -G stringtie_merged.gtf-o ~/ballgown/ERR204916/ERR204916_chrX.gtf ERR204916_chrX.bam

[if !supportLists]8、[endif]用ballgown包進(jìn)行可視化。

通過bioconductor網(wǎng)站安裝需要的R包：

Ballgown、dplyr、DelayedMatrixStats

在R中運行如下程序：

library(ballgown)

setwd("D:/data")#設(shè)定工作目錄

pheno_data

= read.csv("geuvadis_phenodata.csv")#讀取一個儲存了樣本ID和變量的文件，這個文件需要自己建立，打開之后如下圖：

?

pheno_data

bg_chrX

= ballgown(dataDir = "ballgown", samplePattern = "ERR",

pData=pheno_data)#將上述數(shù)據(jù)讀取到ballgown函數(shù)中

bg_chrX

bg_chrX_filt= subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)# texpr函數(shù)是ballgown包中的一個計算轉(zhuǎn)錄本表達(dá)水平的函數(shù)，其默認(rèn)輸出值為FPKM,這一步是將轉(zhuǎn)錄本的表達(dá)水平高于1的篩選出來。

if(!requireNamespace("BiocManager", quietly = TRUE))

install.packages("BiocManager")

BiocManager::install("DelayedMatrixStats")

library(DelayedMatrixStats)#安裝R包

results_transcripts=stattest(bg_chrX_filt,feature="transcript",covariate="sex",adjustvars= c("population"), getFC=TRUE, meas="FPKM")

results_genes

=

stattest(bg_chrX_filt,feature="gene",covariate="sex",adjustvars

= c("population"), getFC=TRUE, meas="FPKM")#這里的stattest函數(shù)是為了保持實驗變量唯一，因為我們的樣本中含有兩個變量：性別和種群，而我們在此處只是想要探究一下性別對于基因表達(dá)的影響，因此需要矯正一下種群對基因和轉(zhuǎn)錄本水平變化的影響。

results_transcripts

= data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt),

results_transcripts)#將基因的名字和ID加到輸出數(shù)據(jù)的表格中，為了下面輸出為csv格式的文件做準(zhǔn)備。

library(dplyr)

results_transcripts= arrange(results_transcripts,pval)

results_genes

= arrange(results_genes,pval)#將輸出結(jié)果按照p-value進(jìn)行排序

write.csv(results_transcripts,"chrX_transcript_results.csv", row.names=FALSE)

write.csv(results_genes,

"chrX_gene_results.csv", row.names=FALSE)#結(jié)果輸出，其中row.names=FALSE是為了去除第一列函數(shù)自帶的編號，這些編號沒有什么實際意義。

subset(results_transcripts,results_transcripts$qval<0.05)

write.csv(subset(results_transcripts,results_transcripts$qval<0.05),"diff_transcripts.csv", row.name=TRUE)

write.csv(subset(results_genes,results_genes$qval<0.05),

"diff_genes.csv", row.name=TRUE)#將qvalue<0.05的轉(zhuǎn)錄本篩選出來并且輸出為csv格式。The q-value is an adjusted p-value,

taking in to account the false discovery rate (FDR)R中有一個函數(shù)可以單獨完成這個轉(zhuǎn)換，stattest函數(shù)的輸出結(jié)果中直接進(jìn)行了這一轉(zhuǎn)換。

[if !supportLists]9、[endif]將上述數(shù)據(jù)做成更加形象的圖片

tropical=c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')

palette(tropical)#定義顏色

fpkm= texpr(bg_chrX,meas="FPKM")

fpkm

= log2(fpkm+1)#將fpkm轉(zhuǎn)化成對數(shù)形式主要是為了能夠在一張圖片中進(jìn)行展示。

boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')#col參數(shù)是為了設(shè)置顏色，las是為了設(shè)置數(shù)據(jù)間隔，ylab是為了設(shè)置軸標(biāo)簽。

運行后得到下圖：

關(guān)于 表示基因表達(dá)水平和FPKM的計算公式:

?

?

?