水稻转录组分析Hisat2+Stringtie+Deseq2

关于流程:

 Nature Communications 文章

Gaining comprehensive biological insight into the tranome by performing a broad-spectrum RNA-seq analysis

文章从流行的39个工具中组合了120种RNA-seq分析方法,并对分析结果的精确度、效率和一致性三个层次进行了评估。

在文章的最后作者给出了ballgown自己推荐的RNAseq 分析方法  分析流程如图所示,回帖用HISAT2,组装和定量用StingTie,差异计算选择DESeq2。



软件:

安装miniconda

wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh

sh Miniconda3-latest-Linux-x86_64.sh

source ~/.bashrc


hisat2 和 stringtie 都可以在anaconda上直接安装,deseq2可以使用新的Bioconductor的安装方式:

conda activate

conda create -n hisat2

conda activate hisat2

conda install hisat2


conda activate

conda create -n  stringtie

conda activate stringtie

conda install stringtie


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

install.packages("BiocManager")

BiocManager::install("DESeq2", version="3.8")

下载wget http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip

#解压unzip Trimmomatic-0.39.zip

conda安装

condainstalltrimmomatic




预处理:

和其他数据处理一样,先把sra文件转为fastq,可以使用seqtk对fastq进行自动的trim,但还是要先用FASTQC先看一下质量,太差的数据最好还是不要用。

fastqc -o outdir -t 8 input.fastq

#-o 输出文件目录

#-t 进程数

1建索引

/share/home/caas_zhangZP/miniconda3/envs/hisat2/bin/extract_splice_sites.py rice.gtf > rice.ss

/share/home/caas_zhangZP/miniconda3/envs/hisat2/bin/extract_exons.py rice.gtf >rice.exon

/share/home/caas_zhangZP/miniconda3/envs/hisat2/bin/hisat2-build -p 50 --exon rice.exon --ss rice.ss rice.fa rice

2比对

hisat2 --dta --p 8 -x rice.fa -1 1C.R1.fastq.gz -2 1C.R2.fastq.gz -S 1C.sam

循环比对

for i in {1..18};

do

hisat2  --dta  --p 8 -x rice.fa -1 R${i}_1.out.fq -2 R${i}_2.out.fq -S  R${i}.sam

done

批量比对

vi hisat2.sh

for i in *1.clean.fastq.gz;

do

i=${i%_1.clean.fastq.gz*};

echo ${i};

nohup hisat2 -p 4 --dta -x rice.fa -1

${i}_1.clean.fastq.gz -2 ${i}_1.clean.fastq.gz -S ${i}.sam &

done

3sam转bam并且排序

samtools sort -@ 12 -o SRR${i}.sort.bam SRR${i}.sam

或者分两步

#将sam文件转化为二进制的bam文件,节约空间和处理时间

samtools view --threads 6 -m 2G -bS 1C.sam > 1C.bam

#将bam文件排序(必须进行这一步):

samtools sort --threads 5 -m 2G 1C.bam -o 1C.sorted.bam

批量处理

vi view.sh

for i in *.sam;

do

i=${i%.sam*};

echo ${i};

nohup samtools view -bhS  -t ~/Rna-seq-SRP200940/GRCh37/genome.fa.fai -@ 2 -o ${i}.bam ${i}.sam &

done

vi sort.sh

for i in *.bam;

do

i=${i%.bam*};

echo ${i};

nohup samtools sort -@ 4 -o ${i}.sort.bam ${i}.bam &

done


6定量

对样本进行组装

比对上的reads将会被呈递给StringTie进行转录本组装,StringTie单独的对每个样本进行组装,在组装的过程中顺带估算每个基因及isoform的表达水平

stringtie -p 8 -G rice.gtf -o file.gtf file.bam

使用脚本批量组装

vi stringtie1.sh

for i in *sort.bam;

do

i=${i%.sort.bam*};

echo ${i};

nohup stringtie -p 4 -G  ~/Rna-seq-SRP200940/GRCh37/genome.gtf

-o ${i}.gtf -l  ${i} ${i}.sort.bam &

done


7.将所有转录本合并

所有的转录本都被呈递给StringTie的merge函数进行merge,这一步是必须的,因为有些样本的转录本可能仅仅被部分reads覆盖,无法被StringTie组装出来。merge步骤可以创建出所有样本里面都有的转录本

创建mergelist

ls /.gtf >mergelist.txt

转录本合并

stringtie --merge -p 8 -G rice.gtf  -o stringtie_merged.gtf  mergelist.txt

检测相对于注释基因组,转录本的组装情况

gffcompare -r rice.gtf -G -o merged stringtie_merged.gtf

重新组装转录本并估算基因表达丰度

stringtie –e –B -p 4 -G stringtie_merged.gtf -o ballgown/SRR9228751/SRR9228751.gtf SRR9228751.sort.bam

批量处理

vi chongzuzhuang.sh

for i in *.sort.bam;

do

i=${i%.sort.bam*};

echo ${i};

nohup stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/${i}/${i}.gtf  ${i}.sort.bam &

done


stringtie -p 8 -G genome.gtf -o file.gtf file.bam

#单个样本转录本计算,生成含有初始表达量的gtf文件

ls /.gtf >mergelist.txt

#将所有样本gft文件位置生成一个list

stringtie --merge -p 8 -G genome.gtf -o stringtie_merged.gtf mergelist.txt

#合并所有样本的转录本,生成完整的的转录本拼接结果gtf文件

stringtie -e -B -p 8 -G stringtie_merged.gtf -o ballgown/file1/file1.gtf file1.bam

#重新组装转录本并估算基因表达丰度,并为ballgown创建读入文件


使用tringtie拼接序列

for i in {1..18};

do

    #samtools sort -@ 8 -o R${i}.bam R${i}.sam

    stringtie  -p 8  --rf -G /8Tdata/liuhailong/genome/Sus_scrofa.Sscrofa11.1.97.gtf  -o /8Tdata/liuhailong/final_data/ballgown/R${i}/R${i}.gtf -B -e  R${i}.bam

done












参考资料:

https://www.jianshu.com/p/216ce84a6220

https://www.jianshu.com/p/cf4c5776fc79

https://cloud.tencent.com/developer/article/2032035

https://cloud.tencent.com/developer/article/2056795

https://www.jianshu.com/p/2cbc739669d5

https://blog.csdn.net/flashan_shensanceng/article/details/125548474

https://blog.csdn.net/weixin_39628160/article/details/110667978

https://www.cnblogs.com/triple-y/p/14246809.html

https://www.jianshu.com/p/b86e5598468b

featureCounts和StringTie是两种广泛使用的生物信息学工具,主要用于RNA-seq数据的分析,但它们的功能和重点有所不同。下面是两者的主要区别:

featureCounts

featureCounts是一个由Aaron Lun等人开发的程序,主要功能是基于比对后的reads(如BAM文件)计数到基因组特征上,如基因、外显子、启动子或基因组的固定bin。它通常用于以下目的:

读段计数:featureCounts 将reads分配给基因组中的已知特征,如基因或转录本,并输出每种特征的读段计数,这对于后续的差异表达分析至关重要。

兼容多种比对器:它可以处理由多种比对器生成的输出,包括STAR、Bowtie、TopHat等。

灵活性:用户可以指定多种参数,如是否考虑多映射reads、是否计数到外显子或基因等。

并行化:featureCounts 支持多线程,可以在多个CPU核心上同时运行,加快处理速度。

StringTie

StringTie是一个用于从RNA-seq数据中组装转录本的工具,其主要功能包括:

转录本组装:StringTie 可以从比对后的reads中重新构建转录本,包括已知的和新发现的转录本。这意味着它可以识别出新的剪接变体或未知转录本。

定量表达:除了组装之外,StringTie 还可以估计每个转录本的表达量,通常使用FPKM(fragments per kilobase of transcript per million mapped reads)作为单位。

更新注释:StringTie 可以利用RNA-seq数据更新现有的基因注释,这对于那些基因注释可能不完整或过时的物种尤为重要。

与其他工具的整合:StringTie 的输出可以被其他工具如Ballgown或Cuffdiff2用于进一步的分析,如差异表达分析。

总结

featureCounts 主要关注的是基于已有注释的特征计数,而 StringTie 更侧重于从头组装转录本和更新基因注释。

如果你的目标是进行差异表达分析且已有高质量的注释集,featureCounts 可能是一个直接的选择。

如果你对发现新的转录本或更新基因注释感兴趣,那么 StringTie 或许更适合你。

两者经常在RNA-seq管道中组合使用,其中StringTie可用于组装和更新注释,而featureCounts则用于基于更新后的注释进行读段计数。


featureCounts 和 StringTie 都可以用于RNA-seq数据分析,但在得到FPKM(Fragments Per Kilobase of transcript per Million mapped reads)方面,它们的处理方式有本质的不同。

featureCounts

featureCounts 主要是用来对RNA-seq数据进行基因或转录本级别的计数,它并不直接计算FPKM或类似标准化表达量。featureCounts 输出的是原始的计数(counts),即有多少reads被映射到了特定的基因或转录本上。然而,这些原始计数可以被进一步处理来计算FPKM或其他标准化表达量指标,例如通过使用额外的软件包如DESeq2或edgeR。这是因为featureCounts只关注于计数,而不涉及基于长度和测序深度的标准化。

StringTie

相比之下,StringTie 不仅可以组装转录本,还可以直接计算每个转录本的FPKM。StringTie 在组装转录本的同时,考虑到转录本的长度和整个文库的测序深度,从而计算出FPKM值。这使得StringTie能够提供标准化后的表达量,可以直接用于比较不同样本之间转录本的表达水平,无需额外的标准化步骤。

总结

featureCounts 提供的是未标准化的读段计数,这些计数可以被其他工具用来计算FPKM或类似标准化表达量。

StringTie 在转录本组装的过程中直接计算FPKM,提供标准化后的表达量。

当使用featureCounts时,通常会采用额外的分析步骤(如使用DESeq2或edgeR)来进行标准化和差异表达分析。而StringTie则提供了一站式的解决方案,包括转录本的组装和表达量的标准化。

在实际应用中,选择哪种工具取决于你的具体需求:如果你需要高度精确的转录本组装和直接的FPKM值,StringTie可能是更好的选择;而如果你只需要基因级别的计数并且计划使用其他工具进行标准化和差异表达分析,那么featureCounts就足够了。

或者分两步

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容