当前位置: 首页 > news >正文

做网站已经不行视频推广方案模板

做网站已经不行,视频推广方案模板,外贸网站如何做推广,兰州官网排名推广在正式分析之前,对于数据的处理是至关重要的,这种重要性是体现在很多方面,其中有一点是要求分析者采用正确的数据类型。 对于芯片数据,原始数据进行log2处理之后可以进行很多常见的分析,比如差异分析、热图、箱线图、…

在正式分析之前,对于数据的处理是至关重要的,这种重要性是体现在很多方面,其中有一点是要求分析者采用正确的数据类型。

对于芯片数据,原始数据进行log2处理之后可以进行很多常见的分析,比如差异分析、热图、箱线图、PCA分析、生存分析、模型构建,聚类分析和相关性分析等。

对于转录组数据,在上述的常见分析中只有差异分析时需要采用count值其他的分析是需要采用log2后的cpm,tpm,fpkm,rpkm数据

除了上述的常规分析,在使用不同R包进行分析之前务必浏览一下输入数据的要求

那么芯片数据还好说,毕竟后续进行log2处理后就可以做很多分析。但是转录组数据的可选项就比较多了。

但目前常用的只有tpm和cpm

count数据转换为cpm数据非常简单

# exprSet是count表达矩阵
# 一句代码搞定
exprSet = log2(edgeR::cpm(exprSet)+1)

比较难的就是count数据转换为tpm数据,因此搬运了常规的流程和R包的方法,做个对比

首先要去获取基因长度文件,因为后续需要用这个数据去矫正基因长度。

网址:https://gdc.cancer.gov/about-data/gdc-data-processing/gdc-reference-files

方法一:常规方法-来自生信技能树

1、获取gencode.v36.annotation并处理

if(!require(GenomicFeatures))BiocManager::install("GenomicFeatures")
library(GenomicFeatures)
txdb <- makeTxDbFromGFF("gencode.v36.annotation.gtf.gz",format="gtf")
# 获取每个基因id的外显子数据
exons.list.per.gene <- exonsBy(txdb,by="gene")
# 对于每个基因,将所有外显子减少成一组非重叠外显子,计算它们的长度(宽度)并求和
exonic.gene.sizes <- sum(width(GenomicRanges::reduce(exons.list.per.gene)))
# 得到geneid和长度数据
gfe <- data.frame(gene_id=names(exonic.gene.sizes),length=exonic.gene.sizes)
head(gfe)[1:5,1:2]
#                               gene_id length
# ENSG00000000003.15 ENSG00000000003.15   4536
# ENSG00000000005.6   ENSG00000000005.6   1476
# ENSG00000000419.13 ENSG00000000419.13   1207
# ENSG00000000457.14 ENSG00000000457.14   6883
# ENSG00000000460.17 ENSG00000000460.17   5970
save(gfe,file = "gfe.Rdata")

2、使用TCGA-HNSC数据

load("hnsc_exp.Rdata")
head(hnsc)[1:3,1:3]
#                      TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15                         2090                         1680                         5433
# ENSG00000000005.6                             0                            0                            0
# ENSG00000000419.13                         2098                         3872                         2240
identical(rownames(hnsc),rownames(gfe))
#[1] TRUE 行名是能够对上的

使用了TCGA-HNSC中的count数据,检测了一下count数据和下载的基因信息的顺序是一致的。

3、曾老师写的代码进行count/tpm转化

load("gfe.Rdata")
#提取基因长度列
effLen = gfe$length
#转化
Counts2TPM <- function(counts, effLen){rate <- log(counts) - log(effLen)denom <- log(sum(exp(rate)))exp(rate - denom + log(1e6))
}
hnsc_tpm_raw <- apply(hnsc, 2, Counts2TPM, effLen = effLen)
head(hnsc_tpm_raw)[1:3,1:3]
#                      TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15                     31.13713                     15.74248                     72.37614
# ENSG00000000005.6                       0.00000                      0.00000                      0.00000
# ENSG00000000419.13                    117.46366                    136.35307                    112.14231

4、加载既往处理好的数据进行对比

load("hnsc_tpm.Rdata")
head(hnsc_tpm)[1:3,1:3]
#                       TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# ENSG00000000003.15                     31.13713                     15.74248                     72.37614
# ENSG00000000005.6                       0.00000                      0.00000                      0.00000
# ENSG00000000419.13                    117.46366                    136.35307                    112.14231

结果是一致的

方法二:R包 DGEobj.utils
load("gfe.Rdata")
#提取上面处理好的基因长度列
effLen = gfe$length
library(DGEobj.utils)
CC_res <- convertCounts(exp,unit = "TPM", # CPM、FPKM、FPK 或 TPMgeneLength = effLen, #这里还是需要下在基因长度数据log = FALSE, #默认为False,设为TRUE时返回Log2值normalize = "none", #默认为none,调用edgeR的calcNormFactors()进行标准化prior.count = NULL #为避免取0的对数,对每个观测值添加平均count。仅当 log = TRUE 时使用,
)
head(CC_res)[1:3,1:3]
#                   TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
#ENSG00000000003.15                     31.13713                     15.74248                     72.37614
#ENSG00000000005.6                       0.00000                      0.00000                      0.00000
#ENSG00000000419.13                    117.46366                    136.35307                    112.14231

与前面得到的分析结果是一致的。

需要下载基因长度数据,但是前期处理完之后后面可以很方便的转化为各种想要的格式(CPM、FPKM、FPK 或 TPM)。

方法三:R包 IOBR
load("hnsc_exp.Rdata")
exp <- hnsc
library(IOBR)
exp_tpm <- count2tpm(exp, #表达矩阵, 行为基因idType = "Ensembl", #"Ensembl" "ENTREZ","SYMBOL" 
)
head(exp_tpm)[1:3,1:3]
#         TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# TSPAN6                     31.89988                     16.16771                     73.85059
# TNMD                        0.00000                      0.00000                      0.00000
# DPM1                       15.65888                     18.22162                     14.88931
# 可以看到这个函数还自动把gene_id给转换成了symbol# 为了更好的对比,我们也把hnsc_tpm_raw中的gene_id转换成symbol
# 使用小洁老师的trans_exp_new函数
library(tinyarray)
a <- trans_exp_new(hnsc_tpm_raw)
head(a)[1:3,1:3]
#           TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1                       0.000000                    0.0000000                     0.000000
# WASH7P                        1.150477                    0.8180029                     2.191641
# MIR6859-1                     0.993794                    0.0000000                     4.443138# 排个序再看一下
p <- identical(rownames(a),rownames(exp_tpm));p
if(!p) {s = intersect(rownames(a),rownames(exp_tpm))a = a[s,]exp_tpm = exp_tpm[s,]
}
head(a)[1:3,1:3]
#          TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1                     0.000000                   0.00000000                     0.000000
# WASH7P                      1.150477                   0.81800290                     2.191641
# FAM138A                     0.000000                   0.03486849                     0.000000head(exp_tpm)[1:3,1:3]
#          TCGA-UF-A7JF-01A-11R-A34R-07 TCGA-CN-4725-01A-01R-1436-07 TCGA-D6-6827-01A-11R-1915-07
# DDX11L1                      0.00000                   0.00000000                      0.00000
# WASH7P                      22.63327                   5.95811827                     15.88102
# FAM138A                      0.00000                   0.03581036                      0.00000

IOBR包中的count2tpm函数只能进行tpm转化(github上搜了这个R包内容似乎没有转化为其他数据格式的函数了)。用默认的方式进行运算得到的结果存在一定的偏差,而且我个人觉得这个偏差有点大... 但是我暂时不知道是什么原因?是内置的基因长度顺序有问题?还是我某个参数设的不对?求高手指点~

综合上述分析,暂时还是选择常规方法DGEobj.utils R包中的convertCounts函数吧~

参考资料:

1、https://hbctraining.github.io/Training-modules/planning_successful_rnaseq/lessons/sample_level_QC.html

2、https://rdrr.io/cran/DGEobj.utils/man/convertCounts.html

3、https://rdrr.io/github/IOBR/IOBR/man/count2tpm.html

4、生信技能树推文:https://mp.weixin.qq.com/s/IUV9dSbRBK1nvetixKOCRw

致谢:感谢曾老师,小洁老师以及生信技能树团队全体成员。

:若对内容有疑惑或者有发现明确错误的朋友,请联系后台(欢迎交流)。更多内容可关注公众号:生信方舟

- END -

http://www.ritt.cn/news/4514.html

相关文章:

  • 成都网站排名优化开发哈尔滨最新信息
  • 免费网站建设软件百度商家入驻怎么做
  • 国外网站代做武汉seo软件
  • 惠州建设工程质量监督站网站安徽网站seo公司
  • 国外做彩票网站客服站长工具手机综合查询
  • 网站建设与运营及营销服务营销型网站建设的步骤流程是什么
  • 建设银行网站用户注册不了今日最新国际新闻头条
  • 网站建设行业资讯广州网络营销产品代理
  • 网页设计教程ui河北seo基础教程
  • 电视台网站如何做新闻报道互联网营销工具有哪些
  • 做网站有必要吗百度竞价推广开户
  • 卖磁铁的网站怎么做win10一键优化工具
  • 成都教育网站建设公司价格网上营销模式
  • 免费网站开发框架seo关键词外包公司
  • 学做视频t的网站windows优化大师官网
  • cms系统都有哪些郑州seo网络推广
  • 网站vip怎么做百度网页链接
  • 微信网站开发系统seo关键词排名优化联系方式
  • 山东苹果网站建设方案网站优化推广方案
  • 网页前端技术企业关键词优化专业公司
  • 上海建站正规电商培训学校排名
  • 网络销售模式 自建网站sem是什么意思
  • 西安网站建设昆奇百度收录查询代码
  • 网站建设门户网站自动秒收录工具
  • 网站建设公司特色电商运营seo
  • 中跃建设集团网站吗唐山百度搜索排名优化
  • 做服装批发的网站哪个比较好在线识别图片
  • 网站文件命名规则郑州百度推广代运营
  • wordpress后台好用优化排名 生客seo
  • 广告传媒公司网站博客营销