参考文献：
Pertea M, Kim D,Pertea G M, et al. Transcript-level expression analysis of RNA-seq experimentswith HISAT, StringTie and Ballgown.[J]. Nature Protocols, 2016, 11(9):1650.
1.clean reads质量检验：
用fastqc检验clean reads质量，代码如下：
fastqc -o 【待测序列所在目录】--extract -f fastq *.fq.gz（序列名称）
质检结果：
所有序列的Kmer Content均为FAIL，Per Base Sequence Content、Per Tile sequence quality、sequence duplication levels至少有一个WARN。
Per Base Sequence Content图中四条线交织，表明存在overrepresented sequence：
猜测：转录组中，表达量高的reads被识别为overrepresented sequence，属于正常现象。
参考：http://www.cnblogs.com/longjianggu/p/5078782.html
2.制作索引（indexes）：
在GENOSCOPE下载油菜基因组序列文件（.fa.gz）和基因组注释文件（.gff3.gz）：http://www.genoscope.cns.fr/brassicanapus/data/
hisat2-build /genome/GCF_000686985.1_Brassica_napus_assembly_v1.0_genomic.fna（基因组序列文件目录） brassica_tran（索引名称及输出目录）
3.序列比对（alignment）：
原始代码：
hisat2 -p（线程数） 16 --dta -x indexes/brassica_tran（索引文件目录） -1 /data/CleanData/$reads1（reads1目录） -2 /data/CleanData/$reads2（reads2目录） -S /hisat2_results/$【sample name】.sam（输出文件名称及目录）
每个样本都需通过上述代码进行比对，封装成BASH脚本运行：
vi alignment（新建空白文档alignment，以下代码在vi中输入）
#! /bin/bash
# run some hisat2 alignments
while read line
do
reads1=${line}_1.fq.gz
reads2=${line}_2.fq.gz
hisat2 -p 16 --dta -x /indexes/brassica_tran -1/data/CleanData/$reads1 -2 /data/CleanData/$reads2 -S/hisat2_results/${line}.sam
done </files/samples.txt（将样本名称按行写在位于/files目录下的samples.txt文件中）
<Esc>:wq保存退出
chmod +x alignment（赋予可执行权限）
nohup ./alignment & 后台运行脚本，输出结果将储存在生成的nohup.log文件中
4.samtools排序及格式转换
原始代码：
samtools sort -@ 16 -o /data/bamfiles/【sample name】.bam（bam文件名称及输出目录） /data/hisat2_results/【sample name】.sam（sam文件名称及输出目录）
封装成BASH脚本：
#! /bin/bash
# samtools sort
while read line
do
samtools sort -@ 16 -o/data/bamfiles/${line}.bam /data/hisat2_results/${line}.sam
done < /files/samples.txt
5.序列初组装（assembly）
原始代码：
stringtie -p 16 -G /genes/brassica.gff（基因组注释文件） -o /data/transcripts/【sample name】.gtf（输出，比对结果，gtf文件） -l 【sample name】（命名规则） /data/bamfiles/【sample name】.bam（输入，bam文件所在目录）
封装成BASH脚本：
#! /bin/bash
# stringtie 1st step
while read line
do
stringtie -p 16 -G/genes/brassica.gff -o /data/transcripts/${line}.gtf -l ${line}/data/bamfiles/${line}.bam
done < /files/samples.txt
6.Merge
stringtie --merge -p 16 -G /genes/brassica.gff -o stringtie_merged.gtf/data/transcripts/mergelist.txt
gffcompare检测组装结果：
gffcompare -r /genes/brassica.gff -G -o merged stringtie_merged.gtf
7.计算表达量并输出成ballgown格式
原始代码：
stringtie -e -B -p 16 -G /data/transcripts/stringtie_merged.gtf（用merge生成的gtf文件代替基因组注释） -o /data/ballgown/$【sample name】/$【sample name】.gtf（输出为ballgown所需的输入格式） /data/bamfiles/【sample name】.bam（输入，bam文件）
封装成BASH脚本：
#! /bin/bash
# stringtie 3rd step
while read line
do
stringtie -e -B -p 16 -G/data/transcripts/stringtie_merged.gtf -o /data/ballgown/${line}/${line}.gtf/data/bamfiles/${line}.bam
done < /files/samples.txt
8.ballgown分析差异表达基因（在R中进行）
>install.packages(‘devtools’)
>source('http://www.bioconductor.org/biocLite.R')
>biocLite(c(’ballgown’, 'genefilter’,'dplyr’,'devtools’))
>library(ballgown)
>library(genefilter)
>library(dplyr)
>library(devtools)
>file <- ballgown(dataDir='/data/ballgown',samplePattern=【sample名前缀】）输入基因表达数据
>samplesNames(file)查看ballgown文件中各样本的排列顺序
>pData(file) <- read.csv(‘/data/geuvadis_phenodata.csv’）（导入自制的表型数据）
>filefilt <-subset(file,'rowVars(texpr(file))>1',genomesubset=TRUE)（过滤掉表达差异较小的基因）
>diff_genes <- stattest(filefilt,feature='gene',covariate=【自变量】,adjustvars=【无关变量】,meas='FPKM')（统计差异表达的基因）
将差异表达基因按pval从小到大排序，并写入csv文件：
>diff_genes <- arrange(diff_genes,pval)
>write.csv(diff_genes,'/data/diff_genes',row.names=FALSE)
————————————————
版权声明：本文为CSDN博主「ygyxl」的原创文章，遵循CC 4.0 BY-SA版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/ygyxl/article/details/74375031
步均做完，不知道第8步.ballgown分析差异表达基因（在R中进行）怎么办了，一直报错。比如>install.packages(‘devtools’)
报错如下
