代码文件备份 | 6-30 组学数据差异分析

于 2022-06-30 19:38:13 发布
阅读量255 收藏 1
点赞数
分类专栏： 生物信息学基础分析 文章标签： 生物信息学
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/weixin_40640700/article/details/125546952
版权 
setwd("E://组学数据分析实战(3.21)")
phe <- read.csv("Example_metaData.csv", as.is=TRUE)
count<-rep(1,dim(phe)[1])
phe<-cbind(phe,count)
library(ggplot2)
p1 <- ggplot(phe, aes(Species,weight=count,fill=Sex)) +
  geom_bar(color = "black", width = .7, position = 'dodge')+
  labs(y="Sample number",x="Species",title = "the number of samples for each spacies")+
  theme(plot.title = element_text(hjust = 0.5))
p1

#数据结构不对#重新看一下不同物种period的数量

human_time<-as.data.frame(table(phe[phe$Species=="Human",]$Predicted.period))
human_da1<-cbind(Species=rep("Human",dim(human_time)[1]),human_time)
macaque_time<-as.data.frame(table(phe[phe$Species=="Macaque",]$Predicted.period))
macaque_dat1<-cbind(Species=rep("Macaque",dim(macaque_time)[1]),macaque_time)
plotData2<-rbind(human_da1,macaque_dat1)
colnames(plotData2)<-c("Species","predicted.period","count")
plotData2$predicted.period<-as.factor(plotData2$predicted.period)

subtypeOrder<-c("2","3","4","5","6","7","8","9","10","11","12","13")
phe$Predicted.period<-droplevels(factor(phe$Predicted.period,levels=subtypeOrder)) #the level is not right

library(ggplot2)
p2 <- ggplot(plotData2, aes(Species, weight = count,fill= predicted.period)) +
  geom_bar(color = "black", width = .7, position = 'stack') +
  labs(y="Sample number",x="Species",title = "the number of samples for each spacies\n in each predicted period")+
  theme(plot.title = element_text(hjust = 0.5))
p2

#####################################################################################
#for all expressed genes
rpkm<-read.table("Example_RPKM.txt")
phe <- read.csv("Example_metaData.csv", as.is=TRUE)
annova_data<-cbind(phe[,c(2,4,11,15)],t(rpkm))
outcome<-as.matrix(annova_data[,4:26105])

annova_data$Species<-as.factor(annova_data$Species)
#model<-manova(cbind(apply(outcome,2,as.numeric))~Species+Sex+Predicted.period+RIN, data = annova_data)
model<-manova(cbind(annova_data$`ENSG00000226376.1|RP11-29B9.1`,annova_data$`ENSG00000216998.1|CYP2AC1P`,annova_data$`ENSG00000175093.4|SPSB4`)~Species+Sex+Predicted.period+RIN, data = annova_data)
summary(model)

#####################################################################################
rpkm<-read.table("Example_RPKM.txt")
phe <- read.csv("Example_metaData.csv", as.is=TRUE)

#人剔除2，3，4，6，11
#猴剔除4，10
table(phe[phe$Species=="Human",]$Predicted.period)
table(phe[phe$Species=="Macaque",]$Predicted.period)
#看一下对应的label有哪些？（这样做有无意义）
phe[phe$Species=="Human"&phe$Predicted.period=="2",]$Sample
phe[phe$Species=="Human"&phe$Predicted.period=="3",]$Sample
phe[phe$Species=="Human"&phe$Predicted.period=="4",]$Sample
phe[phe$Species=="Human"&phe$Predicted.period=="6",]$Sample
phe[phe$Species=="Human"&phe$Predicted.period=="7",]$Sample
phe[phe$Species=="Human"&phe$Predicted.period=="11",]$Sample
phe[phe$Species=="Macaque"&phe$Predicted.period=="4",]$Sample
phe[phe$Species=="Macaque"&phe$Predicted.period=="7",]$Sample
phe[phe$Species=="Macaque"&phe$Predicted.period=="10",]$Sample
filter<-c("HSB112.DFC","HSB148.DFC","HSB113.DFC","HSB150.DFC","HSB153.DFC","HSB103.DFC","HSB114.DFC","HSB149.DFC",
          "HSB107.DFC","HSB159.DFC","HSB92.DFC","HSB98.DFC","HSB194.DFC","HSB175.DFC","RMB209.DFC","RMB208.DFC","RMB200L.DFC","RMB233.DFC",
          "RMB291.DFC","RMB212.DFC","RMB213.DFC","RMB214.DFC","RMB215.DFC")
filter_count<-
filter_data<-rpkm[,!colnames(rpkm)%in%filter]
filter_phe<-phe[!phe$Sample%in%filter,]
dim(filter_phe)
table(filter_phe[filter_phe$Species=="Human",]$Predicted.period)
table(filter_phe[filter_phe$Species=="Macaque",]$Predicted.period)

write.csv(filter_phe,"filter_phe.csv")
write.csv(filter_data,"filter_data.csv")

human_time<-as.data.frame(table(filter_phe[filter_phe$Species=="Human",]$Predicted.period))
human_da1<-cbind(Species=rep("Human",dim(human_time)[1]),human_time)
macaque_time<-as.data.frame(table(filter_phe[filter_phe$Species=="Macaque",]$Predicted.period))
macaque_dat1<-cbind(Species=rep("Macaque",dim(macaque_time)[1]),macaque_time)
plotData2<-rbind(human_da1,macaque_dat1)
colnames(plotData2)<-c("Species","predicted.period","count")
plotData2$predicted.period<-as.factor(plotData2$predicted.period)

#subtypeOrder<-c("2","3","4","5","6","7","8","9","10","11","12","13")
plotData2$Predicted.period<-droplevels(plotData2$Predicted.period) #the level is not right


library(ggplot2)
p2 <- ggplot(plotData2, aes(Species, weight = count,fill= predicted.period)) +
  geom_bar(color = "black", width = .7, position = 'stack') +
  labs(y="Sample number",x="Species",title = "the number of samples for each spacies\n in each predicted period\n(after filter)")+
  theme(plot.title = element_text(hjust = 0.5))
p2

########################################################################################
#pca ananlysis
#首先去除方差为0的行
filter_data2<-filter_data[!apply(filter_data,1,var)%in%c(0,NA),]
pca <- prcomp(t(filter_data2), center=TRUE, scale=TRUE) #注意这里需要转置

pcaPlotDat <- cbind(filter_phe, pca$x[,1:2])
ggplot(pcaPlotDat) + geom_point(aes(x=PC1, y=PC2, shape=Species, color=factor(Predicted.period)), size=3)

#this region maybe can be transferred 
#does the code from the teacher is wrong?
time_gene<-pca$rotation[order(abs(pca$rotation[,"PC1"]),decreasing=TRUE),1:2][1:30,] #正负两个方向变化最小的基因
species_gene<-pca$rotation[order(abs(pca$rotation[,"PC2"]),decreasing=TRUE),1:2][1:30,]
write.csv(rbind(time_gene,species_gene),"pca_specific_gene.csv")

#################################################################

filter_phe<-read.csv("filter_phe.csv")
filter_phe<-filter_phe[,-1]
filter_data<-read.csv("filter_data.csv")
rownames(filter_data)<-filter_data[,1]
filter_data2<-filter_data[,-1]
#plot the interest gene 
scatter_gene<-function(gene){
  plotDat <- cbind(filter_phe[,1:15],gene=t(filter_data2[grep(gene,rownames(filter_data2)),])[,1])
  SpeciesCols <- c("red", "green")
  names(SpeciesCols) <- c("Human", "Macaque")
  SpeciesColors <- SpeciesCols[plotDat$Species]
  SpeciesPchs <- c(20, 17)
  names(SpeciesPchs) <- c("Human", "Macaque")
  SpeciesPointTypes <- SpeciesPchs[plotDat$Species]
  plot(log2(plotDat$Days), plotDat$gene, xlab="Log2(days)", ylab=gene, col=SpeciesColors, pch=SpeciesPointTypes,main = gene)
  legend("topright", inset=0, title="Species", c("Human","Macaque"), pch=SpeciesPchs, col=SpeciesCols)
}
scatter_gene("HLF")
scatter_gene("MEX3A")
scatter_gene("RNF187")
scatter_gene("TERB1")
#wo cuo le!kaolv bazhe sizhangtu huazaiyiqi
pca$rotation[order(pca$rotation[,"PC1"]),1:2][1:30,] 
#似乎看不到比较明显的特征，是我理解错了，还是真的老师写错了。


##################################################################
library(DESeq2)
count <- read.table("Example_count.txt", as.is=TRUE, header=TRUE, sep="\t")
filter<-c("HSB112.DFC","HSB148.DFC","HSB113.DFC","HSB150.DFC","HSB153.DFC","HSB103.DFC","HSB114.DFC","HSB149.DFC",
          "HSB107.DFC","HSB159.DFC","HSB92.DFC","HSB98.DFC","HSB194.DFC","HSB175.DFC","RMB209.DFC","RMB208.DFC","RMB200L.DFC","RMB233.DFC",
          "RMB291.DFC","RMB212.DFC","RMB213.DFC","RMB214.DFC","RMB215.DFC")
filter_count<-count[,!colnames(count)%in%filter]
phe<-read.csv("filter_phe.csv")
phe$Predicted.period <- factor(phe$Predicted.period)
phe$Species <- factor(phe$Species)
dds <- DESeqDataSetFromMatrix(countData = filter_count,
                              colData = phe,
                              design= ~ Species + Predicted.period)

dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, name="Species_Macaque_vs_Human")
sigGenes <- subset(res, padj < 0.01 & abs(log2FoldChange) > 1) #提取到差异表达基因
filter_rpkm<-read.csv("filter_data.csv")
rownames(filter_rpkm)<-filter_rpkm$X
filter_rpkm<-filter_rpkm[,-1]
sigGeneDat <- filter_rpkm[rownames(sigGenes),]
#现在想法就是绘制热图
periodCols <- colorRampPalette(c("yellow", "green", "blue", "purple"))(5)
names(periodCols) <- levels(phe$Predicted.period)
periodColors <- periodCols[as.character(phe$Predicted.period)] #时期
speciesCols<-c("red","blue")
names(speciesCols)<-levels(phe$Species)
speciesColors<-speciesCols[as.character(phe$Species)]
color_label<-cbind(as.character(phe$Species),as.character(phe$Predicted.period))
rownames(color_label)<-colnames(sigGeneDat)
colnames(color_label)<-c("Species","Period")
annotation_col<-as.data.frame(color_label)
library(pheatmap)
pdf("sigGene_heatmap_EuclideanDist5.pdf", 10, 12)
pheatmap(sigGeneDat,scale ="row",annotation_col=annotation_col,color = colorRampPalette(colors = c("blue","white","red"))(100),
       cutree_cols =2 )
dev.off()
#我不明白为什么会报错
#我明白了报错在数据格式上
dim(sigGeneDat)
dim(annotation_col)

pheatmap(sigGeneDat,scale ="row",annotation_col=annotation_col,color = colorRampPalette(colors = c("blue","white","red"))(100),
         cutree_cols =2 )


##########################################
#绘制火山图
res <- results(dds, name="Species_Macaque_vs_Human")
#sigGenes <- subset(res, padj < 0.01 & abs(log2FoldChange) > 1)
volplotDat<-as.data.frame(subset(res,select=c("log2FoldChange","padj")))
volplotDat2<-na.omit(volplotDat)

cut_off_pvalue = 0.01  #统计显著性
cut_off_logFC = 1           #差异倍数值

# 根据阈值参数，上调基因设置为‘up’，下调基因设置为‘Down’，无差异设置为‘Stable’，并保存到change列中
volplotDat2$change = ifelse(volplotDat2$padj< cut_off_pvalue & abs(volplotDat2$log2FoldChange) >= cut_off_logFC, 
                        ifelse(volplotDat2$log2FoldChange> cut_off_logFC ,'Macaque specific','Human specific'),
                        'Stable')
volplotDat2<-cbind(volplotDat2,gene=row.names(volplotDat2))
library(ggplot2)
p3 <- ggplot(
  # 数据、映射、颜色
  volplotDat2, aes(x =log2FoldChange, y = -log10(padj), colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#B22222", "#3CB371","#63B8FF"))+
  # 辅助线
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_pvalue),lty=4,col="black",lwd=0.8) +
  # 坐标轴
  labs(x="log2(fold change)",
       y="-log10 (p.adj)",title = "Volcano Plot(Macaque vs Human)")+
  theme_gray()+
  # 图例
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank())
p3

library(ggrepel)
volplotDat2$label <- ifelse(volplotDat2$padj<cut_off_pvalue& abs(volplotDat2$log2FoldChange) >= 5,
                        as.character(volplotDat2$gene), "")

p3+ geom_label_repel(data = volplotDat2, aes(x = log2FoldChange, 
                                         y = -log10(padj), 
                                         label = label),
                     size = 3, box.padding = unit(0.5, "lines"),
                     point.padding = unit(0.8, "lines"), 
                     segment.color = "black", 
                     show.legend = FALSE)

#scatter_gene("GATS")

par(mfrow=c(2,2))
scatter_gene("BEX5")
scatter_gene("TRIB3")
scatter_gene("HSP90AA5P")
scatter_gene("ZNF33CP")

###################################################################
#analyse first and then illustrate the process
library(DESeq2)
resultsNames(dds)
pe_13_5<- results(dds, name="Predicted.period_13_vs_5")
pe_12_5<- results(dds, name="Predicted.period_12_vs_5")
pe_9_5<- results(dds, name="Predicted.period_9_vs_5")
pe_8_5<- results(dds, name="Predicted.period_8_vs_5")

#接下来一个个分析
sigGenes_5_VS13 <- as.data.frame(subset(pe_13_5, padj < 0.01 & log2FoldChange < -2))
sigGenes_5_VS12 <- as.data.frame(subset(pe_12_5, padj < 0.01 & log2FoldChange < -2))
sigGenes_5_VS9 <- as.data.frame(subset(pe_9_5, padj < 0.01 & log2FoldChange < -2))
sigGenes_5_VS8 <- as.data.frame(subset(pe_8_5, padj < 0.01 & log2FoldChange < -2))
p5_genes <- Reduce(intersect, list(rownames(sigGenes_5_VS12),rownames(sigGenes_5_VS13),rownames(sigGenes_5_VS9),rownames(sigGenes_5_VS8)))

write.csv(p5_genes,"p5_genes.csv")
#sigGenes_5<-rbind(sigGenes_5_VS12,sigGenes_5_VS13,sigGenes_5_VS8,sigGenes_5_VS9)
table(duplicated(rownames(sigGenes_5)))
sigGenes_5<-sigGenes_5[duplicated(rownames(sigGenes_5)),]
#中间的intersect
sigGenes_13 <- as.data.frame(subset(pe_13_5, padj < 0.01 & log2FoldChange >2))
sigGenes_12 <- as.data.frame(subset(pe_12_5, padj < 0.01 & log2FoldChange >2))
sigGenes_9 <- as.data.frame(subset(pe_9_5, padj < 0.01 & log2FoldChange >2))
sigGenes_8 <- as.data.frame(subset(pe_8_5, padj < 0.01 & log2FoldChange >2))


#取他有而其他数据集中没有的,不断的筛
del1<-intersect(rownames(sigGenes_13),rownames(sigGenes_12))
del2<-intersect(rownames(sigGenes_13),rownames(sigGenes_9))
del3<-intersect(rownames(sigGenes_13),rownames(sigGenes_8))
del13<-unique(c(del1,del2,del3))
sig13<-rownames(sigGenes_13)[!row.names(sigGenes_13)%in%del13]
write.csv(sig13,"p13_genes.csv")


del4<-intersect(rownames(sigGenes_12),rownames(sigGenes_9))
del5<-intersect(rownames(sigGenes_12),rownames(sigGenes_8))
del12<-unique(c(del1,del4,del5))
sig12<-rownames(sigGenes_12)[!row.names(sigGenes_12)%in%del12]
write.csv(sig12,"p12_genes.csv")


del6<-intersect(rownames(sigGenes_9),rownames(sigGenes_8))
del9<-unique(c(del2,del4,del5))
sig9<-rownames(sigGenes_9)[!row.names(sigGenes_9)%in%del9]
write.csv(sig9,"p9_genes.csv")


del8<-unique(c(del3,del5,del6))
sig8<-rownames(sigGenes_8)[!row.names(sigGenes_8)%in%del8]
write.csv(rownames(sigGenes_8),"p8_genes.csv")

scatter_gene("SLF1")
scatter_gene("SLC35F3")
scatter_gene("APBA1")
scatter_gene("HGF")
scatter_gene("CLDND1")


siggene_Period<-c(p5_genes,sig8,sig9,sig12,sig13)
filter_rpkm<-read.csv("filter_data.csv")
rownames(filter_rpkm)<-filter_rpkm$X
filter_rpkm<-filter_rpkm[,-1]
sigGeneDat <- filter_rpkm[siggene_Period,]

#现在想法就是绘制热图
phe<-read.csv("filter_phe.csv")
phe<-phe[,-1]
phe$Predicted.period<- factor(phe$Predicted.period,levels = c("5","8","9","12","13"))
phe<-phe[order(phe$Predicted.period),]
periodCols <- colorRampPalette(c("green", "red"))(5)
names(periodCols) <- levels(phe$Predicted.period)
periodColors <- periodCols[as.character(phe$Predicted.period)] #时期
speciesCols<-c("red","blue")
names(speciesCols)<-levels(phe$Species)
speciesColors<-speciesCols[as.character(phe$Species)]
color_label<-cbind(as.character(phe$Species),as.character(phe$Predicted.period))
#color_label<-cbind(as.character(phe$Predicted.period))
rownames(color_label)<-colnames(sigGeneDat)
colnames(color_label)<-c("Species","Period")
#colnames(color_label)<-c("Period")
annotation_col<-as.data.frame(color_label)
annotation_row<-as.data.frame(c(rep("period 5",length(p5_genes)),
                                    rep("period 8",length(sig8)),
                                    rep("period 9",length(sig9)),
                                    rep("period 12",length(sig12)),
                                    rep("period 13",length(sig13))))
colnames(annotation_row)<-"Period_Genes"
rownames(annotation_row)<-rownames(sigGeneDat)
newcols<-c("RMB202.DFC","HSB154.DFC","HSB178.DFC","HSB96.DFC","HSB97.DFC",  
           "RMB193L.DFC","RMB228.DFC","RMB207.DFC","HSB139.DFC",
           "HSB122.DFC","RMB201L.DFC","RMB216.DFC","HSB119.DFC","HSB124.DFC","HSB127.DFC","RMB188L.DFC",
           "RMB199L.DFC","RMB219.DFC","HSB123.DFC","HSB126.DFC","HSB130.DFC","HSB135.DFC","HSB136.DFC", 
           "HSB145.DFC","RMB160.DFC","RMB161.DFC","RMB177R.DFC","RMB196.DFC" )
sigGeneDat<-sigGeneDat[,newcols]


ann_colors = list(
  Period = c(`5`="#CD1076",`8`="#ADFF2F",`9`="#FF8C00",`12`="#228B22",`13`="#1E90FF"),
  Species= c(Human="#DDA0DD",Macaque="#FFFF00"),
  Period_Genes = c(`period 5`="#CD1076",`period 8`="#ADFF2F",`period 9`="#FF8C00",`period 12`="#228B22",`period 13`="#1E90FF")
)
head(ann_colors)

#group_df = data.frame(Groups=as.factor(rep(c("Control", "Treated"), c(5,5))))
#ann_colors = list(
#  Groups = c(Control="black", Treated="white"))

library(pheatmap)
#pheatmap(sigGeneDat,scale ="row",annotation_col=annotation_col,color = colorRampPalette(colors = c("blue","white","red"))(100),
#         cutree_cols =2 )

#这个效果并不是我想要的,有没有其他改进的方法
pheatmap(sigGeneDat,scale ="row",show_rownames = F,show_colnames=F,cluster_rows=F,cluster_cols=F,annotation_row = annotation_row,annotation_col=annotation_col,color = colorRampPalette(colors = c("#006400","white","#FF1493"))(100),
         annotation_colors = ann_colors )

########################################################################
library(plyr)
library(AnnotationHub)	
library(org.Hs.eg.db) 
library(clusterProfiler)
library(dplyr)
library(ggplot2)
library(gridExtra)

res <- results(dds, name="Species_Macaque_vs_Human")
sigGenes_humam <- rownames(as.data.frame(subset(res, padj < 0.01 & log2FoldChange < -2)))
sigGenes_macaque <-rownames(as.data.frame(subset(res, padj < 0.01 & log2FoldChange > 2)))

sig_human_change<-read.table("human_symbol.csv")
#f1<-c(sigGenes_macaque) #要求输入的ensemble id
#f1<-c(sigGenes_humam)
#library(stringr)
#f<-as.data.frame(str_split_fixed(f1,"\\|",2))
f1=sig_human_change
#EG2Ensembl=toTable(org.Hs.egSYMBOL)	 #将ENTREZID和ENSEMBL对应的数据存入该变量
EG2Ensembl=toTable(org.Hs.egENSEMBL)
#f=f1$V2
f=f1$V1
geneLists=data.frame(ensembl_id=f)
results=merge(geneLists,EG2Ensembl,by='ensembl_id',all.x=T)
id=na.omit(results$gene_id)  #提取出非NA的ENTREZID

#par(mfrow=c(2,3))
#GO富集分析
ego <- enrichGO(OrgDb="org.Hs.eg.db", gene = id, ont = "MF", pvalueCutoff = 0.01, readable= TRUE) 
p1<-dotplot(ego,showCategory=10,title="Enrichment GO Top10(MF)") 
p2<-barplot(ego, showCategory=20,title="EnrichmentGO(MF)",cex.axis=0.5) 

#plotGOgraph(ego)
ego2<-enrichGO(OrgDb="org.Hs.eg.db", gene = id, ont = "CC", pvalueCutoff = 0.01, readable= TRUE) 
p3<-dotplot(ego2,showCategory=10,title="Enrichment GO Top10(CC)") 
p4<-barplot(ego2, showCategory=20,title="EnrichmentGO(CC)")

ego3<-enrichGO(OrgDb="org.Hs.eg.db", gene = id, ont = "BP", pvalueCutoff = 0.05, readable= TRUE) 
p5<-dotplot(ego3,showCategory=10,title="Enrichment GO Top10(BP)")+scale_y_discrete(labels = function(x) str_wrap(x, width = 50))
p6<-barplot(ego3, showCategory=20,title="EnrichmentGO(BP) for LINE")

# grid.arrange(p1,p2,p3,p4,p5,p6,ncol=3, nrow = 2)
# grid.arrange(p1,p2,ncol=2, nrow = 1)
# grid.arrange(p3,p4,ncol=2, nrow = 1)
# grid.arrange(p5,p6,ncol=2, nrow = 1)

#KEGG分析#
ekk <- enrichKEGG(gene= id,organism  = 'hsa', qvalueCutoff = 0.05)	 #KEGG富集分析
dotplot(ekk,font.size=8,title="KEGG pathway Enrichment")	# 画气泡图


enrichplot<-function(genes,a){
  f1<-c(genes) #要求输入的ensemble id
  library(stringr)
  f<-as.data.frame(str_split_fixed(f1,"\\|",2))
  EG2Ensembl=toTable(org.Hs.egSYMBOL)	 #将ENTREZID和ENSEMBL对应的数据存入该变量
  f=f$V2
  geneLists=data.frame(symbol=f)
  results=merge(geneLists,EG2Ensembl,by='symbol',all.x=T)
  id=na.omit(results$gene_id)  #提取出非NA的ENTREZID
  ego3<-enrichGO(OrgDb="org.Hs.eg.db", gene = id, ont = "BP", pvalueCutoff = 0.05, readable= TRUE) 
  filename<-paste(a,"plot.jpg",sep = "_")
  p<-dotplot(ego3,showCategory=10,title="Enrichment GO Top10(BP)")+scale_y_discrete(labels = function(x) str_wrap(x, width = 50))
  jpeg(file = filename)
  print(p)
  dev.off()
  }

par(mfrow=c(2,2))
enrichplot(sigGenes_humam,"human_enrich")
#p1<-dotplot(pA,showCategory=10,title="Enrichment GO Top10(BP)")+scale_y_discrete(labels = function(x) str_wrap(x, width = 50))

enrichplot(sigGenes_macaque,"macaque_enrich")
enrichplot(p5_genes,"period_young")
enrichplot(c(sig8,sig9,sig12,sig13),"period_old")

#########################################################################
#vennplot
install.packages("UpSetR")
library(UpSetR)
listInput <- list(sigGenes_humam,sigGenes_macaque,p5_genes,sig8,sig9,sig12,sig13)
names(listInput)<-c("human","macaque","period 5","period 8","period 9","period 12","period 13")
upset(fromList(listInput), order.by = "freq",nsets = 7,nintersects = NA)
intersect(sigGenes_macaque,p5_genes)
intersect(sigGenes_humam,p5_genes)
intersect(sigGenes_macaque,sig13)
intersect(sigGenes_humam,sig13)
par(mfrow=c(2,2),cex=0.5)
#scatter_gene("FNDC11")
scatter_gene("FKBP7")
#scatter_gene("LRRC3")
#scatter_gene("CKAP2L")
#scatter_gene("NEUROG2")
#scatter_gene("LL0XNC01-116E7.1")
#scatter_gene("KB-51A8.1")
#scatter_gene("SDC1")
scatter_gene("RP11-247C2.2")
#scatter_gene("TMEM95")
#scatter_gene("NOX5")
scatter_gene("AVP")
scatter_gene("TMPRSS3")
今天也是个妖精头子呀
关注
  • 0
    点赞
  • 踩
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值