本文参考文献https://www.nature.com/articles/nprot.2016.095#materials进行bulk RNA-seq数据分析
首先进行原始数据的下载，当然因为文章中已经有了链接进行数据的下载，因此直接使用wget命令就能完成所需数据的下载，但是也要学会基本的测序数据、参考基因组、注释文件、index索引文件的下载，当然也可以自己建立index文件。
raw data下载（找了2个小鼠的测序结果）：
fastq-dump --split-3 --gzip SRR3589959
fastq-dump --split-3 --gzip SRR3589960
加上--split-3之后, 会把原来双端拆分成两个文件,但是原来单端并不会保存成两个文件. 还有你用--gzip就能输出gz格式, 能够节省空间的同时也不会给后续比对软件造成压力,比对软件都支持，就是时间要多一点。
标准基因组的数据下载：
注释文件的下载：
Index的下载，自己的电脑没法建立index有现成的index可以去https://ccb.jhu.edu/software/hisat2/index.shtml下载：
关于HISAT2官网上的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
当然也可以自己建立索引文件：
Hisat2-build reference genome
Index文件不会用不就狗带了吗，注意index文件的格式，其实就是只有电脑能看你不能看的二进制，所以千万不要尝试打开，否则会是一堆乱码，感觉自己弄错了emmm，打开index文件是一些.ht2格式的文件，参考基因组有几条染色体就有几个这样的文件，但是运行HISAT2的时候千万不要在文件名后面加上后缀！要不然会报错的~
分析步骤：
[if !supportLists]1、[endif]首先是对测序数据质量的检测——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、去除接头和低质量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]进行第二次质检：
fastqc-o ~/rnaseq/data/secondqc -t 6 ERR188044_chrX_1_val_1.fq.gzERR188044_chrX_2_val_2.fq.gz
结果和原来差不多，因此便用了原来的reads进行比对。
[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
以上参数都可以通过hisat2 –help查到。
4、将sam文件转换成bam文件并进行排序，至于为什么要进行排序是因为打开其中一个sam文件就可以看到人家已经指明数据是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]进行转录本的组装。
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]将组装后的转录本merge一下：
stringtie --merge -p8 -G chrX.gtf -o stringtie_merged.gtf ~/rnaseq/chrX_data/mergelist.txt
其中mergelist.txt就是上一步得到的文件名字组成的纯文本格式的文件，当然也可以将上述文件名直接打上去，以空格隔开。下图为mergelist.txt文件打开的样子：
[if !supportLists]6、[endif]将merge后的注释文件与标准基因组的注释文件进行对比（这一步可以不做）：
gffcompare -rchrX.gtf -G -o merged stringtie_merged.gtf
[if !supportLists]7、[endif]再以merge后的注释文件作为参考进行转录本的组装并输出为ballgown可以识别的格式（因为在一个实验中，我们往往有很多样本，而不同样本的转录本往往是不同的，而我们需要让他们在近似相同的情况下进行组装，因此再用merge后的注释文件作为参考进行组装）：
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包进行可视化。
通过bioconductor网站安装需要的R包：
Ballgown、dplyr、DelayedMatrixStats
在R中运行如下程序：
library(ballgown)
setwd("D:/data")#设定工作目录
pheno_data
= read.csv("geuvadis_phenodata.csv")#读取一个储存了样本ID和变量的文件，这个文件需要自己建立，打开之后如下图：
pheno_data
bg_chrX
= ballgown(dataDir = "ballgown", samplePattern = "ERR",
pData=pheno_data)#将上述数据读取到ballgown函数中
bg_chrX
bg_chrX_filt= subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)# texpr函数是ballgown包中的一个计算转录本表达水平的函数，其默认输出值为FPKM,这一步是将转录本的表达水平高于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函数是为了保持实验变量唯一，因为我们的样本中含有两个变量：性别和种群，而我们在此处只是想要探究一下性别对于基因表达的影响，因此需要矫正一下种群对基因和转录本水平变化的影响。
results_transcripts
= data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt),
results_transcripts)#将基因的名字和ID加到输出数据的表格中，为了下面输出为csv格式的文件做准备。
library(dplyr)
results_transcripts= arrange(results_transcripts,pval)
results_genes
= arrange(results_genes,pval)#将输出结果按照p-value进行排序
write.csv(results_transcripts,"chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes,
"chrX_gene_results.csv", row.names=FALSE)#结果输出，其中row.names=FALSE是为了去除第一列函数自带的编号，这些编号没有什么实际意义。
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的转录本筛选出来并且输出为csv格式。The q-value is an adjusted p-value,
taking in to account the false discovery rate (FDR)R中有一个函数可以单独完成这个转换，stattest函数的输出结果中直接进行了这一转换。
[if !supportLists]9、[endif]将上述数据做成更加形象的图片
tropical=c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)#定义颜色
fpkm= texpr(bg_chrX,meas="FPKM")
fpkm
= log2(fpkm+1)#将fpkm转化成对数形式主要是为了能够在一张图片中进行展示。
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')#col参数是为了设置颜色，las是为了设置数据间隔，ylab是为了设置轴标签。
运行后得到下图：
关于 表示基因表达水平和FPKM的计算公式:
