关于流程:
Nature Communications 文章
文章从流行的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就足够了。
或者分两步