论文
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结尾，这个信息没啥用，给去掉
代码语言：javascript复制new.names <- colnames(my.counts)
new.names <- gsub(new.names,pattern = ".eff.counts", replacement = "")
colnames(my.counts) <- new.names
把第一列转录本id作为行名，再把数据去个整
代码语言：javascript复制rownames(my.counts) <- my.counts[,1]
my.counts.round <- round(my.counts[,-1])
对基因进行过滤，标准是最小的count数是10，转录本最少在5个样本中的count数都大于10
代码语言：javascript复制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
代码语言：javascript复制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")
这个是初步处理好了，后续还需要其他处理敬请期待
image.png
示例数据和代码等整套流程都写完了再分享
