学习Nature正刊论文中eQTL分析前对基因表达量的预处理

小明的数据分析笔记本

论文

Graph pangenome captures missing heritability and empowers tomato breeding

https://www.nature.com/articles/s41586-022-04808-9

论文中的方法部分写到

RNA-seq data were quantified as transcripts per million(TPM). Genes with TPM values of > 0.5 were defined as expressed. Only genes expressed in at least 100 accessions were retained for the downstream analyses. Raw expression levels were normalized with quantile-quantile normalization. To remove potential batch effects and confounding factors affecting gene expression, the probabilistic estimation residuals method were applied with the top four factors as covariates.

表达量数据可以在论文中提供的链接处下载

读取数据，这里表达量文件还挺大的，可以用data.table这个R包的fread函数读取

library(data.table)
library(tidyverse)

dat<-fread("D:/Jupyter/practice/WGCNA_nature_tomato/expression_raw_332.txt")
dim(dat)
dat[1:6,1:4]

这里的数据行是样本，列是基因，首先做个转置

dat %>% 
  column_to_rownames("target_id") %>% 
  t() %>% 
  as.data.frame() -> dat.t

第一步

依据表达量对数据进行过滤，至少在100个样本中TPM值大于0.5

dat.t[rowSums(dat.t > 0.5) >= 100,] -> dat.t.filter
dim(dat.t.filter)

第二步

标准化

preprocessCore::normalize.quantiles(dat.t.filter %>% as.matrix()) %>% 
  as.data.frame() -> dat.t.filter.norm

colnames(dat.t.filter.norm)<-colnames(dat.t.filter)
rownames(dat.t.filter.norm)<-rownames(dat.t.filter)

dat.t.filter.norm[1:5,1:5]

dat.t.filter.norm %>% 
  rownames_to_column("geneid") %>% 
  write_tsv("D:/Jupyter/practice/WGCNA_nature_tomato/dat.t.filter.norm.tsv")

第三步

运行peer

这个需要用到R 3.5版本的环境

代码内容

library(peer)

args <- commandArgs(trailingOnly = TRUE)

expr<-read.table(args[1],row.names=1,header=TRUE,sep="\t")
model = PEER()
PEER_setPhenoMean(model,as.matrix(t(expr)))
PEER_setNk(model,20)
PEER_getNk(model)
PEER_update(model)
factors = as.data.frame(t(PEER_getX(model)))
colnames(factors)<-colnames(expr)
rownames(factors)<-paste0("peer",1:20)
write.table(factors,args[2],quote=F,row.names=T,sep="\t",col.names=TRUE)

time Rscript run_peer.R dat.t.filter.norm.tsv exp.peer.covar.tsv

这一步需要的时间比较长，运行了167分钟

论文中用top4 factors作为协变量，应该是选择输出结果的前5行就可以了

欢迎大家关注我的公众号 小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！

微信扫一扫分享文章