软件官网：
Hisat2： Manual | HISAT2
StringTie：StringTie
文章：Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown | Nature Protocols
建议看保姆级教程：
1. RNA-seq : Hisat2+Stringtie+DESeq2 – 恒诺新知
2. RNA-seq用hisat2、stringtie、DESeq2分析 - 简书
基本用法：
1. 构建参考基因组索引
# 提取剪接位点和外显子信息
extract_splice_sites.py Mus_musculus.GRCm39.104.gtf > Mus_musculus.ss
extract_exons.py Mus_musculus.GRCm39.104.gtf > Mus_musculus.exon
# 建立索引
# 最后的 Mus_musculus.GRCm39_tran 为索引文件前缀
hisat2-build --ss Mus_musculus.ss --exon Mus_musculus.exon Mus_musculus.GRCm39.dna.primary_assembly.fa \
Mus_musculus.GRCm39_tran
# 时间超长，大于12h，建议晚上跑
2. 参考基因组比对
# -x跟索引名前缀，-1，-2跟双端测序文件，-U跟单端测序文件，-S输出为sam格式的文件,-p线程数量
# 我们直接输出为排序好的bam文件
# --dta输出为转录本组装的reads，--summary-file输出比对信息
hisat2 -p 10 --dta -x path/to/Mus_musculus.GRCm39_tran
--summary-file test1_summary.txt
-1 1.fastq-data/test1_R1_rep1.fq.gz
-2 1.fastq-data/test1_R2_rep1.fq.gz
-S test1.sam
3. samtools 对输出 sam 文件排序并转为 bam 文件
# -@为samtools的线程数
samtools sort -@ 10 -o test1.sorted.bam test.sam
4. 转录本组装
# 组装转录本，-p为线程数，-G为组装参考注释文件，-l为输出文件名前缀
# 单个样本运行
stringtie -p 10 -G Mus_musculus.GRCm38.102.gtf
-l test1
-o test1.gtf
test1.sorted.bam
5. 注释文件合并
# 创建 mergelist.txt 文件，指明组装后注释文件的路径
path/to/test1.gtf
path/to/test2.gtf
path/to/test3.gtf
# 合并gtf文件
$ stringtie --merge -p 10 -G ./Mus_musculus.GRCm38.102.gtf
-o stringtie_merged.gtf
mergelist.txt
6. 利用生成的注释文件对转录本进行定量
# 创建一个新的 test1 文件夹，转录本定量结果保存到文件夹中
mkdir test1/
stringtie -p 10 -e -G ./stringtie_merged.gtf
-o test1/test1.gtf
-A test1/gene_abundances.tsv
test1.sorted.bam
# 相应文件夹下生成样本名.gtf和gene_abundances.tsv的两个文件，对应每个样本的 count 值定量结果，我们需要合并到一个文件里。
7. 提取基因定量结果
prepDE.py 需要一个 sample_list，第一列为样本名，第二列为 gtf 文件路径
# sample_list.txt 文件内容如下
test1 path/to/test1/test1.gtf
test2 path/to/test1/test2.gtf
test3 path/to/test1/test3.gtf
test4 path/to/test1/test4.gtf
# 提取合并count结果，-i为输入sample_list
prepDE.py -i sample_list.txt
# 生成gene_count_matrix.csv和transcript_count_matrix.csv文件
8. 选做：提取 FPKM/TPM 或 coverage 结果
需要用到stringtie_expression_matrix.pl，下载地址如下：
rnaseq_tutorial/stringtie_expression_matrix.pl at master · griffithlab/rnaseq_tutorial · GitHub
# 提取TPM
$ ./stringtie_expression_matrix.pl --expression_metric=TPM
--result_dirs='test1_rep1,test1_rep2,test2_rep1,test2_rep2'
--transcript_matrix_file=transcript_tpms_all_samples.tsv
--gene_matrix_file=gene_tpms_all_samples.tsv
# 提取FPKM
./stringtie_expression_matrix.pl --expression_metric=FPKM
--result_dirs='test1_rep1,test1_rep2,test2_rep1,test2_rep2'
--transcript_matrix_file=transcript_fpkms_all_samples.tsv
--gene_matrix_file=gene_fpkms_all_samples.tsv
# 提取coverage
./stringtie_expression_matrix.pl --expression_metric=coverage
--result_dirs='test1_rep1,test1_rep2,test2_rep1,test2_rep2'
--transcript_matrix_file=transcript_coverage_all_samples.tsv
--gene_matrix_file=gene_coverage_all_samples.tsv
# 在当前目录就会生成相应的基因和转录本的tpm、fpkm、coverage 结果
9. DESeq2 差异分析
# 安装DESeq2包
BiocManager::install('DESeq2')
# 加载包
library(DESeq2)
# 设置工作路径
setwd('D:rnaseq')
# 读入counts矩阵
gene_count_matrix <- read.csv("D:/rnaseq/gene_count_matrix.csv",row.names = 1)
count <- gene_count_matrix[rowSums(gene_count_matrix)>0,]
# 构建表型矩阵
colData <- data.frame(row.names = colnames(count),
condition = factor(c(rep('control',2),rep('treat',2)),
levels=c('control','treat'))
)
# 查看
colData
# condition
# test1_rep1 control
# test1_rep2 control
# test2_rep1 treat
# test2_rep2 treat
dds <- DESeqDataSetFromMatrix(countData = count, colData = colData,design = ~ condition)
dds <- DESeq(dds)
res <- results(dds)
diff_res <- as.data.frame(res)
diff_res$gene_name <- rownames(diff_res)
# 输出差异结果
write.table(diff_res,file = 'DESeq2_diff_results.csv',quote = F,sep = ',',row.names = F,col.names = T)
本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)
