1.读取raw count表达矩阵和分组信息

rawcount <- read.table("data/raw_counts.txt",row.names = 1, 
                       sep = "\t", header = T)
group <- read.table("data/filereport_read_run_PRJNA229998_tsv.txt",
                    header = T,sep = "\t", quote = "\"")

提取表达矩阵对应的样本表型信息

library(stringr)
group <- group[match(colnames(rawcount), group$run_accession),
               c("run_accession","sample_title")]
group$sample_title <- str_split_fixed(group$sample_title,pattern  = "_", n=2)[,2]
group
write.table(group,file = "data/group.txt",row.names = F,sep = "\t",quote = F)

这里主要是因为raw count表达矩阵和分组信息并非一一对应，因此需要用match在分组信息里面挑选

2.表达矩阵预处理

keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)##
filter_count <- rawcount[keep,]
filter_count[1:4,1:4]
dim(filter_count)

这里的预处理主要是处理raw count表达矩阵的含有0的行，代码的意思是如果一行的表达量大于0的个数有0.75倍的行数就保留这行，反之就删除这行。

计算counts per millio(cpm) 表达矩阵，这里用到的是edgeR包的cpm函数

library(edgeR)
express_cpm <- cpm(filter_count)
express_cpm[1:6,1:6]
# 保存表达矩阵和分组结果
save(filter_count, express_cpm, group, file = "data/Step01-airwayData.Rdata")