跟着Genes|Genomes|Genetics学数据分析：WGCNA分析前期的数据预处理01

用户7010445

发布于 2023-01-06 20:15:15

5330

发布于 2023-01-06 20:15:15

文章被收录于专栏：小明的数据分析笔记本

论文

Sex-Specific Co-expression Networks and Sex-Biased Gene Expression in the Salmonid Brook Charr Salvelinus fontinalis

数据代码公开

https://github.com/bensutherland/sfon_wgcna

还有wgcna的代码，论文里对方法和结果部分介绍的还挺详细，可以对照着论文然后学习WGCNA的代码

今天的推文我们学习一下wgcna数据分析前对表达量矩阵的的一些预处理代码

表达量矩阵是原始的count值，原始的count值为什么有小数，我查了一下，找到了一个解答

https://hbctraining.github.io/Intro-to-rnaseq-hpc-O2/lessons/05_counting_reads.html

image.png

这里表达量矩阵文件有98M，使用R语言自带的read.csv()函数读取相对挺慢的，可以使用R包readr中的read_csv()函数读取

my.counts <- read_csv(file = "out.matrix.csv")
dim(my.counts)
my.counts[1:6,1:2]

数据集格式行是基因，列是样本

每个样本的名字都以.eff.counts结尾，这个信息没啥用，给去掉

new.names <- colnames(my.counts)
new.names <- gsub(new.names,pattern = ".eff.counts", replacement = "")
colnames(my.counts) <- new.names

把第一列转录本id作为行名，再把数据去个整

rownames(my.counts) <- my.counts[,1]
my.counts.round <- round(my.counts[,-1])

对基因进行过滤，标准是最小的count数是10，转录本最少在5个样本中的count数都大于10

library(edgeR)
my.counts <- DGEList(counts = my.counts.round)
min.reads.mapping.per.transcript <- 10
cpm.filt <- min.reads.mapping.per.transcript / min(my.counts$samples$lib.size) * 1000000
min.ind <- 5

keep <- rowSums(cpm(my.counts)>cpm.filt) >= min.ind
table(keep)
my.counts <- my.counts[keep, , keep.lib.sizes=FALSE]

把数据集转换成cpm的值，用到的是edgeR这个R包

还可以把cpm值进行log2转化

cpm This unit is known as counts per million reads mapped (CPM)

详细可以参考 http://luisvalesilva.com/datasimple/rna-seq_units.html

dim(my.counts)
my.counts <- calcNormFactors(my.counts, method = c("TMM"))
my.counts <- estimateDisp(my.counts)

normalized.output <- cpm(my.counts, normalized.lib.sizes = TRUE, log= F)

write.csv(normalized.output, file = "03_normalized_data/normalized_output_matrix.csv")
normalized.output.log2 <- cpm(my.counts, normalized.lib.sizes = TRUE, log= T, prior.count = 1)
write.csv(normalized.output, file = "03_normalized_data/normalized_output_matrix_log2.csv")

这个是初步处理好了，后续还需要其他处理敬请期待