转录组上游分析,Count计算

本期教程原文链接:转录组定量,最简单的操作,你会吗?

本期教程

第六章 转录本定量分析

定量软件有RSEM,eXpress,salmoe,kallistofeatureCounts。在网络中吗,都有比较详细的教程,大家可以自己去学习。

本教程中推荐使用两种方式获得转录本的表达量,stringtie -eBfeatureCounts,stringtie自带的脚本程序prepDE.py

6.1 Stringtie -eB

Stringtie -eB是通过stringtie组装后的merge.gtf注释信息二次与.bam文件进行转录本表达量的比对,获得转录本的FPKM,此后使用Ballgown 包结合使用,进行后续的分析。详情可看Transcript-level expression analysis of RNA-seq experiments with HISAT, StringTie and Ballgown, 2016, Nat Protoc一文。

此步骤操作后,每个样本会获得新的gtf文件。

stringtie –e –B -p 8 -G stringtie_merged.gtf -o ballgown/ERR188245/ERR188245_chrX.gtf ERR188245_chrX.bam

使用R语言中的ballgown进行分析

pheno_data = read.csv("geuvadis_phenodata.csv")
bg_chrX = ballgown(dataDir = "ballgown", samplePattern = "ERR", pData=pheno_data)
bg_chrX_filt = subset(bg_chrX,"rowVars(texpr(bg_chrX)) >1",genomesubset=TRUE)
results_transcripts = stattest(bg_chrX_filt, feature="transcript",covariate="sex",adjustvars = c("population"), getFC=TRUE, meas="FPKM")
##
# Identify genes that show statistically significant differences between groups. For this we can run the same  function that we used to identify differentially expressed transcripts, but here we set feature="gene" in the stattest command:
results_genes = stattest(bg_chrX_filt, feature="gene", covariate="sex", adjustvars = c("population"), getFC=TRUE, meas="FPKM")

#  Add gene names and gene IDs to the results_transcripts data frame:
results_transcripts = data.frame(geneNames=ballgown::geneNames(bg_chrX_filt), geneIDs=ballgown::geneIDs(bg_chrX_filt), results_transcripts)

#Sort the results from the smallest P value to the largest:
results_transcripts = arrange(results_transcripts,pval)
results_genes = arrange(results_genes,pval)

# Write the results to a csv file that can be shared and distributed:
write.csv(results_transcripts, "chrX_transcript_results.csv", row.names=FALSE)
write.csv(results_genes, "chrX_gene_results.csv", row.names=FALSE) 
# Identify transcripts and genes with a q value <0.05:
subset(results_transcripts,results_transcripts$qval<0.05)
subset(results_genes,results_genes$qval<0.05)

# Make the plots pretty. This step is optional, but if you do run it you will get the plots in the nice colors that we used to generate our figures:

tropical= c('darkorange', 'dodgerblue', 'hotpink', 'limegreen', 'yellow')
palette(tropical)

##  Show the distribution of gene abundances (measured as FPKM values) across samples, colored by sex (Fig. 3).  
fpkm = texpr(bg_chrX,meas="FPKM")
fpkm = log2(fpkm+1)
boxplot(fpkm,col=as.numeric(pheno_data$sex),las=2,ylab='log2(FPKM+1)')

# Make plots of individual transcripts across samples. For example, here we show how to create a plot for the 12th transcript in the data set (Fig. 4). The first two commands below show the name of the transcript (NM_012227) and the name of the gene that contains it (GTP binding protein 6, GTPBP6):
ballgown::transcriptNames(bg_chrX)[12]
##      12 
## "NM_012227"
ballgown::geneNames(bg_chrX)[12]
##      12 
## "GTPBP6"
plot(fpkm[12,] ~ pheno_data$sex, border=c(1,2), main=paste(ballgown::geneNames(bg_chrX)[12],' :', ballgown::transcriptNames(bg_chrX)[12]),pch=19, xlab="Sex", ylab='log2(FPKM+1)')
points(fpkm[12,] ~ jitter(as.numeric(pheno_data$sex)), col=as.numeric(pheno_data$sex))

6.2 Stringtie中的prepDE.py程序

prepDE.py是stringtie软件中自带的获得转录本丰度的脚本,6.1中的HISAT2+Stringtie+Ballgown是一组黄金组合,但是还有很多的局限性。因此,个人建议仍是获得count值,后再进一步的分析,这样的方法更利于下游分析。

具体详情可查看官方文档:http://ccb.jhu.edu/software/stringtie/index.shtml?t=manual

注意

  1. prepDE.py脚本需要输入strintie -e输出的gtf文件
  2. 此外,需要获得gtf文件的路径,sample_lst.txt

本期教程原文链接:转录组定量,最简单的操作,你会吗?

prepDE.py -h
Usage: prepDE.py [options]

Generates two CSV files containing the count matrices for genes and
transcripts, using the coverage values found in the output of `stringtie -e`

Options:
  -h, --help            show this help message and exit
  -i INPUT, --input=INPUT, --in=INPUT
                        a folder containing all sample sub-directories, or a
                        text file with sample ID and path to its GTF file on
                        each line [default: ./]
  -g G                  where to output the gene count matrix [default:
                        gene_count_matrix.csv
  -t T                  where to output the transcript count matrix [default:
                        transcript_count_matrix.csv]
  -l LENGTH, --length=LENGTH
                        the average read length [default: 75]
  -p PATTERN, --pattern=PATTERN
                        a regular expression that selects the sample
                        subdirectories
  -c, --cluster         whether to cluster genes that overlap with different
                        gene IDs, ignoring ones with geneID pattern (see
                        below)
  -s STRING, --string=STRING
                        if a different prefix is used for geneIDs assigned by
                        StringTie [default: MSTRG]
  -k KEY, --key=KEY     if clustering, what prefix to use for geneIDs assigned
                        by this script [default: prepG]
  -v                    enable verbose processing
  --legend=LEGEND       if clustering, where to output the legend file mapping
                        transcripts to assigned geneIDs [default: legend.csv]
  • 运行:
python prepDE.py -i sample_lst.txt
## 或服务器中的Python版本是Python2或python3
prepDE.py -i sample_lst.txt
#or
prepDE.py3 -i sample_lst.txt
  • sample_lst.txt ,根据自己stringtie -e -B的输出路径进行设置
04_Result/Stringtie_eB/SRR6929571/SRR6929571.gtf
04_Result/Stringtie_eB/SRR6929572/SRR6929572.gtf
04_Result/Stringtie_eB/SRR6929573/SRR6929573.gtf
04_Result/Stringtie_eB/SRR6929574/SRR6929574.gtf
04_Result/Stringtie_eB/SRR6929577/SRR6929577.gtf
04_Result/Stringtie_eB/SRR6929578/SRR6929578.gtf
  • 输出结果:

最终获得gene_count_matrix.csv转录本文件,利用此文件即可进行下有分析。

gene_count_matrix.csv
transcript_count_matrix.csv

每个基因间使用,号隔开

## gene_count_matrix.csv
gene_id,SRR6929571,SRR6929572,SRR6929573,SRR6929574,SRR6929577,SRR6929578
gene:Solyc02g160570.1,0,0,0,0,0,0
gene:Solyc02g161280.1,0,0,0,0,0,0
gene:Solyc01g017370.2,0,0,0,0,0,0
MSTRG.28366,3998,4147,4277,3955,4164,2001

## transcript_count_matrix.csv
transcript_id,SRR6929571,SRR6929572,SRR6929573,SRR6929574,SRR6929577,SRR6929578
mRNA:Solyc06g071220.1.1,19,26,41,62,110,33
mRNA:Solyc07g161730.1.1,0,0,0,0,0,0
MSTRG.23978.1,3,0,0,0,0,0
mRNA:Solyc06g062940.5.1,27523,17069,18415,14892,15523,36301
MSTRG.28358.2,13,18,57,32,6,9
mRNA:Solyc05g055535.1.1,0,0,0,0,0,0
MSTRG.28358.1,123,20,410,192,52,134

为方便下一步的分析,可以,号转变成\t分隔符,常说的Tab分隔符。

sed 's/,/\t/g' gene_count_matrix.csv > 01.gene_count.csv

本期教程原文链接:转录组定量,最简单的操作,你会吗?

6.3 featureCounts的使用

featureCounts是subread中脚本,可以使用subread流程进行定量,在这里直接使用前面mapped的bam文件进行转录本定量。

安装:

conda安装

conda install -y subread

源码安装:

官网:[https://sourceforge.net/projects/subread/

wget https://sourceforge.net/projects/subread/files/subread-2.0.3/subread-2.0.3-Linux-x86_64.tar.gz
tar -zxvf subread-2.0.3-Linux-x86_64.tar.gz
cd subread-2.0.3-Linux-x86_64
cd bin/
#
echo 'PATH=$PATH:~/software/subread-2.0.3-Linux-x86_64/bin' >> ~/.bachrc

运行:

使用featureCounts进行定量的bam文件,我建议使用前期使用hisat2。bowtie2,bwa,或是STAR等mapped的文件。

  • featureCounts可以对转录本(trancript_id)进行定量
  • featureCounts也可以对基因(gene_id)进行定量
Usage: featureCounts [options] -a <annotation_file> -o <output_file> input_file1 [input_file2] ... 

featureCounts -T 5 -p -t exon -g transcript_id -a annotation.gtf -o  counts.txt   *.bam

参数:

-T 
	运行线程数
-t 
	设置feature-type,在GTF注释中指定特征类型。如果提供多个类型,它们之间应以','隔开,中间没有空格。默认情况下是 "exon"
-g 
	GTF注释文件需要计算的基因或转录本的表达水平。默认是:gene_id,可选择transcript_id
-a
	提供的注释文件,可以是参考基因组的annotation.gtf,也可以是组转后的gtf文件
-J
	对可变剪切进行计数
-G < string >
	当-J设置的时候,通过-G提供一个比对的时候使用的参考基因组文件,辅助寻找可变剪切
-o < string >
	输出文件的名字,输出文件的内容为read 的统计数目
-O
	允许多重比对,即当一个read比对到多个feature或多个metafeature的时候,这条read会被统计多次
-d < int >
	最短的fragment,默认是50
-D < int >
	最长的fragmen,默认是600

-p
	能用在paired-end的情况中,会统计fragment而不统计read
-B
	在-p选择的条件下,只有两端read都比对上的fragment才会被统计
-C
	在-p选择的条件下,只有两端read都比对上的fragment才会被统计

运行:

featureCounts -T 20 -p -t exon -g transcript_id -a stringtie_merged.gtf -o All.transcript.count.txt *.sort.bam

获得结果:

All.transcript.count.txt
All.transcript.count.txt.summary

本期教程原文链接:转录组定量,最简单的操作,你会吗?

  • count结果

6.4 Htseq定量

HTseq-count也是一个比较常用的软件,与featureCount功能一样用来计数(count)

安装:

conda install -y htseq

htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m intersection-nonempty yourfile_name.bam ~/reference/hisat2_reference/Homo_sapiens.GRCh38.86.chr_patch_hapl_scaff.gtf > counts.txt
usage: htseq-count [options] alignment_file gff_file

-f {sam,bam}  (default: sam)
	#设置输入文件的格式,该参数的值可以是sam或bam。
-r {pos,name}  (default: name)
	#设置sam或bam文件的排序方式,该参数的值可以是name或pos。
	#前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。
	#若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,
	#两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,
	#程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。
	#因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。
	#而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,
	#也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。
-s {yes,no,reverse}  (default: yes) 
    #数据是否来源于链特异性测序,链特异性是指在建库测序时,只测mRNA反转录出的cDNA序列,而不测该cDNA序列反向互补的另一条DNA序列;换句话说就是,链特异性能更准确反映出mRNA的序列信息
-a MINAQUAL (default: 10)
    #忽略比对质量低于此值的比对结果
-t FEATURETYPE 
    #feature type (3rd column in GFF file) to be used, all features of other type are ignored (default, suitable for Ensembl GTF files: exon)
    #程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。
-i IDATTR
    #GFF attribute to be used as feature ID (default, suitable for Ensembl GTF files: gene_id)
    #设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。
-m {union,intersection-strict,intersection-nonempty} (default: union)
	#设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。
-o | --samout
	#输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。
-n 
	#指定多线程,默认是1

运行Htseq-count:

htseq-count -f bam -r name -s no -a 10 -t exon -i transcript_id -m intersection-strict ../../03_MappedFile/Hisat2_Mapped/*.bam ../../02_Geneome_index/ITAG4.1_gene_models.gtf > htseq_counts.txt

6.5 其他流程

salmon 流程

软件官网:https://combine-lab.github.io/salmon/

6.6 Count to FPKM

  • count to FPKM

使用Perl脚本进行转换

## 需要信息
1. 基因名
2. 基因长度
3. count值

使用cut提取信息:

cat All.transcript.count.txt | cut -f 1,6-13 > 01.all.count.txt

运行perl脚本:

perl CountToFPKM.pl 01.all.count.txt > 02.all.FPKM.txt

CountToFPKM.pl脚本:

本期教程原文链接:转录组定量,最简单的操作,你会吗?

往期教程部分内容











往期部分文章

1. 复现SCI文章系列专栏

2. 《生信知识库订阅须知》,同步更新,易于搜索与管理。

3. 最全WGCNA教程(替换数据即可出全部结果与图形)

  • WGCNA分析 | 全流程分析代码 | 代码一

  • WGCNA分析 | 全流程分析代码 | 代码二

  • WGCNA分析 | 全流程代码分享 | 代码三

  • WGCNA分析 | 全流程分析代码 | 代码四

  • WGCNA分析 | 全流程分析代码 | 代码五(最新版本)


4. 精美图形绘制教程

  • 精美图形绘制教程

5. 转录组分析教程

  • 转录组上游分析教程[零基础]

  • 一个转录组上游分析流程 | Hisat2-Stringtie

6. 转录组下游分析

  • 批量做差异分析及图形绘制 | 基于DESeq2差异分析

  • GO和KEGG富集分析

  • 单基因GSEA富集分析

  • 全基因集GSEA富集分析

小杜的生信筆記 ,主要发表或收录生物信息学的教程,以及基于R的分析和可视化(包括数据分析,图形绘制等);分享感兴趣的文献和学习资料!!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.mfbz.cn/a/553714.html

如若内容造成侵权/违法违规/事实不符,请联系我们进行投诉反馈qq邮箱809451989@qq.com,一经查实,立即删除!

相关文章

最新UI发卡盗U,支持多语言,更新UI界面,支持多个主流钱包,附带系统搭建教程

环境&#xff1a;Linux系统 进入宝塔安装环境&#xff1a;Nginx 1.22.1 MySQL 8.0 php7.4 phpMyAdmin 5.2 按照说明去安装环境&#xff0c;如果没有找到MySQL8.0版本去"软件商店"搜索Mysql切换至8.0 1.上传开源源码 2.上传数据库文件 3.上传猴导入数据库文件 4.修…

SpringBoot+FreeMaker

目录 1.FreeMarker说明2.SpringBootFreeMarker快速搭建Pom文件application.properties文件Controller文件目录结构 3.FreeMarker数据类型3.1.布尔类型3.2.数值类型3.3.字符串类型3.4.日期类型3.5.空值类型3.6.sequence类型3.7.hash类型 4.FreeMarker指令assign自定义变量指令if…

HTML学习笔记:(一)基础方法

Html格式 里面文件使用平台为&#xff1a;w3school 1、基础功能&#xff1a; <html><head> <title>这是我的第一个html页面,会显示在浏览器的标题栏中</title> </head> <!--修改背景颜色 --> <body bgcolor"yellow"> …

QtQuick 学习笔记(二)按钮组件

1. QPushButton 功能 按压按钮&#xff0c;用于接受用户点击事件&#xff0c;可显示设定的字符串提示信息&#xff0c;但需要父组件作为容器&#xff0c;多用以执行命令或触发事件 常用函数 QPushButton&#xff1a;&#xff1a;QPushButton&#xff08;const QString &…

RHCE1

unit1.定时任务和延迟任务项目 1.在系统中设定延迟任务要求如下: 在系统中建立easylee用户&#xff0c;设定其密码为easylee 延迟任务由root用户建立 要求在5小时后备份系统中的用户信息文件到/backup中确保延迟任务是使用非交互模式建立 再使用chmod修改权限&#xff1a; 确保…

HarmonyOS真机调试页面运行卡顿/黑屏解决方法,亲测有效

项目场景&#xff1a; 提示&#xff1a;这里简述项目相关背景&#xff1a; 用mate40等发行时间相对较早但系统是HarmonyOS4.0的真机调试 问题描述 提示&#xff1a;这里描述项目中遇到的问题&#xff1a; 程序点击容易卡顿或黑屏 原因分析&#xff1a; CPU兼容问题导致屏幕…

Text2sql的一些技巧

最近看到了一篇关于text2sql的文章&#xff0c;以及一些论文。对使用模型做text2sql给了一些不错的建议。 参考文章&#xff1a;24年大模型潜力方向&#xff1a;大浪淘沙后的Text-to-SQL和Agent - 知乎 论文&#xff1a;https://arxiv.org/pdf/2403.09732.pdf 关于模型的建议 …

STM32H7定时器TIM1-TIM17中断、PWM实现

STM32H7定时器TIM1-TIM17中断、PWM实现 高级定时器硬件框图定时器模式时基输出PWM定时器输入捕获 TIM1-TIM17的中断配置TIM1-TIM17的PWM输出 STM32H7 支持的定时器有点多&#xff0c;要简单的区分下。STM32H7 支持 TIM1-TIM8&#xff0c;TIM12-TIM17 共14 个定时器&#xff0c;…

使用API有效率地管理Dynadot域名,锁定账户中的域名

关于Dynadot Dynadot是通过ICANN认证的域名注册商&#xff0c;自2002年成立以来&#xff0c;服务于全球108个国家和地区的客户&#xff0c;为数以万计的客户提供简洁&#xff0c;优惠&#xff0c;安全的域名注册以及管理服务。 Dynadot平台操作教程索引&#xff08;包括域名邮…

Zynq7000系列中的时钟管理

PS&#xff08;处理系统&#xff09;时钟子系统生成的所有时钟都源自三个可编程PLL&#xff08;锁相环&#xff09;中的一个&#xff1a;CPU、DDR和I/O。时钟子系统的主要组件如图25-1所示。 在正常工作期间&#xff0c;PLL被启用&#xff0c;并由PS_CLK时钟引脚驱动。在启用P…

6.6Python之集合的基本语法和特性

集合&#xff08;Set&#xff09;是Python中的一种无序、不重复的数据结构。集合是由一组元素组成的&#xff0c;这些元素必须是不可变数据类型&#xff0c;但在集合中每个元素都是唯一的&#xff0c;即集合中不存在重复的元素。 集合的基本语法&#xff1a; 1、元素值必须是…

24卫生高级职称报名时间汇总⏰报名全流程

⏰卫生高级职称&#xff08;网上报名&#xff09;时间汇总 ✔️陕西&#xff1a;4月23日-5月24日 ✔️上海&#xff1a;4月23日-5月24日 ✔️重庆&#xff1a;4月23日—5月24日 ✔️黑龙江&#xff1a;4月23日-5月24日 ✔️浙江&#xff1a;4月23日-5月24日 ✔️云南&#xff1…

面试自救指南:女生如何巧答私密问题

在面试过程中&#xff0c;女性应聘者可能会遇到一些私人问题&#xff0c;这些问题可能涉及婚姻、家庭、生育等方面。面对这些问题&#xff0c;如何回答才能既保持真实又不失礼节呢&#xff1f; 当遇到关于婚姻状况的问题时&#xff0c;您可以选择回答&#xff1a;“我目前的婚姻…

【Python深度学习系列】网格搜索神经网络超参数:权重初始化方法(案例+源码)

这是我的第262篇原创文章。 一、引言 在深度学习中&#xff0c;超参数是指在训练模型时需要手动设置的参数&#xff0c;它们通常不能通过训练数据自动学习得到。超参数的选择对于模型的性能至关重要&#xff0c;因此在进行深度学习实验时&#xff0c;超参数调优通常是一个重要的…

2024全新快递平台系统独立版小程序源码|带cps推广营销流量主+前端

本文来自&#xff1a;2024全新快递平台系统独立版小程序源码|带cps推广营销流量主前端 - 源码1688​​​​​ 应用介绍 快递代发快递代寄寄件小程序可以对接易达云洋一级总代快递小程序&#xff0c;接入云洋/易达物流接口&#xff0c;支持选择快递公司&#xff0c;三通一达&am…

【leetcode面试经典150题】57. 环形链表(C++)

【leetcode面试经典150题】专栏系列将为准备暑期实习生以及秋招的同学们提高在面试时的经典面试算法题的思路和想法。本专栏将以一题多解和精简算法思路为主&#xff0c;题解使用C语言。&#xff08;若有使用其他语言的同学也可了解题解思路&#xff0c;本质上语法内容一致&…

电动车违停智能监测摄像机

电动车的普及带来了便利&#xff0c;但也衍生了一些问题&#xff0c;其中最常见的之一就是电动车的违停。电动车的违停不仅会影响交通秩序&#xff0c;还可能对周围环境和行人安全造成影响。为了监测和管理电动车的违停情况&#xff0c;可以使用电动车违停智能监测摄像机。这种…

退市危机袭来,环保行业能否逆境崛起?|中联环保圈

近年来&#xff0c;环保行业风波持续不断&#xff0c;众多环保大公司风险频出。博天环境的退市危机令人感慨&#xff0c;深圳星源因涉嫌信息披露违法违规而被警告退市&#xff0c;更是引发业界震动。 最近三年&#xff0c;证监会办理的上市公司信息披露违法案件多达 397 件&…

Linux内核之virt_to_page实现与用法实例(五十)

简介&#xff1a; CSDN博客专家&#xff0c;专注Android/Linux系统&#xff0c;分享多mic语音方案、音视频、编解码等技术&#xff0c;与大家一起成长&#xff01; 优质专栏&#xff1a;Audio工程师进阶系列【原创干货持续更新中……】&#x1f680; 优质专栏&#xff1a;多媒…

Python 使用 pip 安装 matplotlib 模块(精华版)

pip 安装 matplotlib 模块 1.使用pip安装matplotlib(五步实现):2.使用下载的matplotlib画图: 1.使用pip安装matplotlib(五步实现): 长话短说&#xff1a;本人下载 matplotlib 花了大概三个半小时屡屡碰壁&#xff0c;险些暴走。为了不让新来的小伙伴走我的弯路&#xff0c;特意…
最新文章