HISAT2 - StringTie - DESeq2 pipeline 进行bulk RNA-seq

软件官网：

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     treatdds <- 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)

HISAT2 - StringTie - DESeq2 pipeline 进行bulk RNA-seq相关推荐

HISAT2+StringTie+Ballgown安装及使用流程
HISAT2+StringTie+Ballgown安装及使用流程 2015年Nature Methods上面发表了一款快速比对工具hisat,作为接替tophat和bowtie的比对工具,它具有更快的 ...
转录组分析---Hisat2+StringTie+Ballgown使用
转录组分析---Hisat2+StringTie+Ballgown使用 (2016-10-10 08:14:45) 转载▼ 标签: 生物信息学转录组 1.Hisat2建立基因组索引: First ...
一文掌握RNA seq，RNA seq课程大汇总
RNA测序(RNA-seq)在过往十年里逐渐成为全转录组水平分析差异基因表达和研究mRNA差异剪接必不可少的工具.RNA-seq帮助大家对RNA生物学的理解会越来越全面:从转录本在何时何地转录到RNA ...
HISAT2+STRINGTIE+BALLGOWN 分析转录组数据
师兄推荐这篇文章,按照里面的命令,先做一套转录组分析. 参考文献: Pertea M, Kim D,Pertea G M, et al. Transcript-level expression ana ...
带有RESTEasy + JAXB + Jettison的JSON示例
RESTEasy使用Jettison JSON库在JSON之间来回映射JAXB注释对象. 在本教程中,我们向您展示如何将带有JAXB注释的对象转换为JSON格式并将其返回给客户端. 杰克逊(Jac ...
string deseq2
hisat2+stringtie+deseq2分析RNA-SEQ数据 - 简书 hisat2+stringtie+deseq2 RNAseq(二) - 简书 RNA-seq : Hisat2+Stri ...
重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述）
原文链接: https://www.embopress.org/doi/10.15252/msb.20188746 主编评语这篇文章最好的地方不只在于推荐了工具,提供了一套分析流程,更在于详细介绍了 ...
对一篇单细胞RNA综述的评述：细胞和基因质控参数的选择
原文链接: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6072887 摘要单细胞RNA测序技术的发展加深了我们对于细胞作为功能单元的理解,不仅能基于成 ...
复现原文（一）：Single-cell RNA sequencing of human kidney（step by step）
https://www.nature.com/articles/s41597-019-0351-8 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这 ...

HISAT2 - StringTie - DESeq2 pipeline 进行bulk RNA-seq

HISAT2 - StringTie - DESeq2 pipeline 进行bulk RNA-seq相关推荐

最新文章

热门文章