前天同事帮领导写基金，看文献找思路，突然问我能不能帮她分析一下数据，我看了一下数据是FPKM的文件。第一反应是可以。我就试试了一下。
1.数据下载
a=read.table("GSE115422_BM_EC_9GY_fpkms.txt.gz",sep='\t',quote = "",fill = T, comment.char = "!",header = T) # 提取表达矩阵
#去重复
b <- a[a$gene_name,]
c <- unique(b)
rownames(c) <- c[,1]
c <- c[,-1]
#选择所需数据
a <- c[,-()]
c <- c[,10:15]
dat1<-as.data.frame(c)
dim(dat1)
dat1[1:4,1:4]
#很多表达量为0的样本，直接选择在某个基因如果在3个样品中的表达量为零，则直接舍去。
apply(dat1,1,function(x){sum(floor(x)==0)>3})
dat1<-dat1[!apply(dat1,1,function(x){sum(floor(x)==0)>3}),]
dim(dat1)
head(dat1)
write.csv(dat1,"dat1.csv")
head(dat1)
boxplot(dat1)#差异很大取log归一化
dat3 <- log(dat1)#下游分析的结果有缺失值，故选择log(dat2 + 1)
dat4 <- log(dat1 + 1)
boxplot(dat3)
boxplot(dat4)
write.csv(dat3,"dat3.csv")
write.csv(dat3,"dat4.csv")
2、差异基因分析
library(limma)
group <- c(rep("EC",3),rep("GY_7d",3), rep("GY_1d",3))
head(group)
View(group)
group <- factor(group)
design <- model.matrix(~0 + group)
colnames(design) <- levels(group)
design
contrast.matrix <- makeContrasts(GY_7d - EC,
GY_1d - EC,
GY_7d - GY_1d,
levels=design)
contrast.matrix
#fit <- lmFit(dat3, design)
fit <- lmFit(dat4,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
allDiff1=topTable(fit2,adjust='fdr',coef=1,number=Inf)
allDiff2=topTable(fit2,adjust='fdr',coef=2,number=Inf)
allDiff3=topTable(fit2,adjust='fdr',coef=3,number=Inf)
save(dat4,group,allDiff1,allDiff2,allDiff3,file = "GSE115422_exprset.Rdata")
write.csv(allDiff1, file = "9GY_7d_vs_EC.csv" );
write.csv(allDiff2, file = "9GY_1d_vs_EC.csv" );
write.csv(allDiff3, file = "9GY_7d_vs_9GY_1d.csv" );
3.如何在差异基因里面显示每一个样本的平均表达量？（同事的需求）
我的答案是取交集，差异基因和原来的表达矩阵取交集。
可能管道更简单些。
first <- read.csv("9GY_7d_vs_EC.csv")
second <- read.csv("9GY_1d_vs_EC.csv")
third <- read.csv("9GY_7d_vs_9GY_1d.csv")
four <- read.csv("dat3.csv")
a <- merge(first, four, by="Gene_name",all=F)
b <- merge(second, four,by="Gene_name",all=F)
c <- merge(third, four,by="Gene_name",all=F)
write.csv(a,file = "9GY_7d_vs_EC1.csv")
write.csv(b,file = "9GY_1d_vs_EC1.csv")
write.csv(c,file = "9GY_7d_vs_9GY_1d1.csv")
后边的热图和火山图我也绘制了，GO和KEGG分析也做了一下，他们想要的是Excel可以读取数据，摸索了很久才用很笨的方法做出来，后边再贴出来。
当时用FPKM做差异分析查了一下，多数是可以做。同事要的比较急就按照上边的流程做了一下。
心里不踏实，因为记得生信技能树上的推文多次提到FPKM，RPKM，TPM意义和区别；翻了一下生信技能树问题完美解决，
- TPM值就是RPKM的百分比
- What the FPKM? A review of RNA-Seq expression units
-
Question: Differential expression analysis starting from TPM data
FPKM转换为TPM：
exprSet <- a
fpkmToTpm <- function(fpkm)
{
exp(log(fpkm) - log(sum(fpkm)) + log(1e6))
}
tpms <- apply(exprSet,2,fpkmToTpm)
tpms[1:3,]
colSums(tpms)
下游就是重复分析一遍流程。
