公式介绍

下方公式是FPKM的综述文章写的，参数并不好理解，但是对数量级表示很明白

https://haroldpimentel.wordpress.com/2014/05/08/what-the-fpkm-a-review-rna-seq-expression-units/

要求出FPKM值，需要获得三个参数。
cDNA Fragments：可以理解为单个基因比对到基因组上的reads数，在测序数据里就是count值。HTseq处理后可以直接得到结果。
Mapped Fragments：指每个样品中所有基因比对到基因组上的reads数。也就是用求和函数sum()将单一样品的count全部加起来。注意，由于单位是百万，所以求和后需要除以10^6。
Transcript Length：也就是exon length，是指reads比对到基因外显子上的长度。这个需要找到参考基因组才能获得数据。

R语言脚本

#！/usr/bin/R#根据参考基因组获得每条reads的exon length /home/XXX/RNA-seq_ref/04_36_homo_2019//ref/gencode/gencode.v29.annotation.gtf
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("/home/XXX/RNA-seq_ref/04_36_homo_2019//ref/gencode/gencode.v29.annotation.gtf",format="gtf")
exons_gene <- exonsBy(txdb, by = "gene")
exons_gene_lens <- lapply(exons_gene,function(x){sum(width(reduce(x)))})
n=t(as.data.frame(exons_gene_lens))

这里是引用网上的一段代码，大概意思就是说使用GenomicFeatures这个R包，可以获得exon length。

需要注意的是GenomicFeatures1.7最新版本只支持R 4.0版本，不支持3.6，请大家重新安装好对应版本再安装工具包。

#导入所有reads的count值，header = T保证后续计算不会将列名算入而造成错误。如果file和脚本不在同一个路径，记得添加file的绝对路径。
count <- read.table(file="all_count.txt",header = T)#从第二列开始，将每个样本的count循环转化成FPKM
i <- 2repeat{mapped_reads <- sum(count[1:58721,i])#计算每个样本的mapped reads数,由于count文本最后几行是综述数据，需要舍去，因此选择58721行之前的基因数据。FPKM <- count[1:58721,i]/(10^-9*mapped_reads*n[1:58721,1])#计算FPKM值，10^-9是单位转换后的换算，不能放在分母外的原因是mapped_reads*n[1:58721,1]的乘积很大会报错，因此需要先缩小倍数再进行计算count = data.frame(count[1:58721,],FPKM)#添加到count表格后面ii <- i + 1if(i >41){break}#样品总共有40个样品，我们需要获取从第2列到第41列count值
}#生成表格列名称
count_colname <- read.table("all_count.txt",header = F,nrow = 1,as.is=TRUE)
FPKM_colname <- paste(count_colname[1,],"_FPKM",sep="")
colname <- c(count_colname,FPKM_colname)
col_name <- colname[-which(colname=="gene_FPKM")]#删除gene_FPKM
names(count) <- col_name#生成表格
write.csv(count,"ALL_FPKM.csv",row.names = FALSE)#row.names是为了去除第一列自动生成的行名，同理col.names=FALSE可以去除第一行自动生成的列名
head(read.csv("ALL_FPKM.csv"))

欢迎扫码或加vx：bbplayer2021进群交流

2020.08.14【RNA-Seq流程】丨将HTseq生成的基因COUNT值转换为FPKM值相关推荐

2021.04.22【RNA-seq流程】丨count值转换为FPKM值优化2.0
优化内容解决每次转换需要设置样本数和基因数目实现基因count值与length精准匹配摘要大概半年前,我写过一篇将HTseq生成的基因COUNT值转换为FPKM值文章,用于对count的入门级 ...
2020.08.14日常总结——Trie树的实际应用
Trie树\color{green}{\texttt{Trie树}}Trie树 Trie树是一种非常高效的字符串处理工具,它把字符串全部放在了一棵树上维护,而且通过一定的手段合并树上的节点,这使得它 ...
2020.09.30【RNA-seq流程】丨转录组生信分析全流程
RNA-Seq生信分析全流程摘要第一部分 step.1 下载数据 step.2 数据质控第二部分 step.3序列比对 step.4 计算基因表达量 step.5 插入片段长度检验 step.6 ...
（十三：2020.08.28）CVPR 2015 追踪之论文纲要（译）
CVPR 2020 追踪之论文纲要(修正于2020.08.27) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
（十一：2020.08.28）CVPR 2017 追踪之论文纲要（译）
CVPR 2017 追踪之论文纲要(修正于2020.08.28) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
（九：2020.08.27）CVPR 2019 追踪之论文纲要（译）
CVPR 2019 追踪之论文纲要(修正于2020.08.28) 讲在前面论文目录讲在前面论坛很多博客都对论文做了总结和分类,但就医学领域而言,对这些论文的筛选信息显然需要更加精细的把控,所以自 ...
心音数据库_小V云端数据库 | 2020.9.14—2020.9.18
桂花的芬芳在雨后空气中弥散开来似为湿润的情绪赠予了一丝甜蜜小V云端数据库 2020.9.14-2020.9.18 资讯情报关键词健康.示范.安全 V宝体检,助力成长 2020年9月14日上午 ...
一文掌握RNA seq，RNA seq课程大汇总
RNA测序(RNA-seq)在过往十年里逐渐成为全转录组水平分析差异基因表达和研究mRNA差异剪接必不可少的工具.RNA-seq帮助大家对RNA生物学的理解会越来越全面:从转录本在何时何地转录到RNA ...
PYTHON学习笔记之（一）2020.08
PYTHON学习笔记之(一)2020.08 Python基础数据类型常见的列表.字典,以及元组.集合. 1 列表 list 1.1 列表转换字符串 stu = ['王一', '李二', '张三'] ...

2020.08.14【RNA-Seq流程】丨将HTseq生成的基因COUNT值转换为FPKM值

count转换FPKM介绍

公式介绍

R语言脚本

2020.08.14【RNA-Seq流程】丨将HTseq生成的基因COUNT值转换为FPKM值相关推荐

最新文章

热门文章