Month: April 2018

  • ggplot2绘制富集分析柱状图

    library(“ggplot2”) pathway=read.table(“COS3_kegg_total.xls”,header=T,sep=”\t”) pathbar = ggplot(pathway,aes(x=Pathway,y=-1*log10(Pvalue))) pathbar + geom_bar(stat=”identity”) pathbar + geom_bar(stat=”identity”) + coord_flip() pathbar+geom_bar(stat=”identity”,color=”red”,fill=”blue”)+coord_flip() pathbar+geom_bar(stat=”identity”,aes(fill=-1*log10(Pvalue)))+coord_flip() porder=factor(as.integer(rownames(pathway)),labels=pathway$Pathway) pathbar+geom_bar(stat=”identity”,aes(x=porder,fill=-1*log10(Pvalue)))+coord_flip() pathbar+geom_bar(stat=”identity”,aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip() porder=factor(as.integer(rownames(pathway)),labels=rev(pathway$Pathway)) pathbar+geom_bar(stat=”identity”,aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip() pqbar=pathbar+geom_bar(stat=”identity”,aes(x=rev(porder),fill=-1*log10(Pvalue)))+coord_flip()+theme(legend.position=”none”)+labs(title=”Top20 of pathway enrichment”,y=expression(-log[10](Pvalue))) pqbar pqbar+geom_hline(yintercept=2,color=c(“red”),linetype=4)

  • pie图绘制

    x

  • awk命令将fastq文件转换成fasta

    awk ‘{if(NR%4==1||NR%4==2){print $0}}’ test.fq | sed ‘s/^@/>/g’ > myfile.fasta 有时候,比如说我们需要做一个clustalw,我们需要将多行的fasta文件(multiple line fasta)格式文件,转换成单行的fasta文件(single line fasta)序列文件,怎么办呢。 巧用awk + printf一行搞定。 awk ‘/^>/ { print n $0;} !/^>/ {printf “%s”, $0, n=”\n”} END {print “”}’ test.fa  

  • WGCNA(my project)

    setwd(“D:/Ma_transcription/WGCNA/”) database <- read.table(file = “COS_deseq_counts_normalized.txt”, sep = “\t”, header = T, row.names = 1, stringsAsFactors = F) library(reshape2) library(‘WGCNA’) enableWGCNAThreads()#打开多线程 WGCNA_matrix = t(database[order(apply(database,1,mad), decreasing = T)[1:10000],]) subname=sapply(colnames(database),function(x) strsplit(x,”_”)[[1]][1]) datTraits = data.frame(gsm=names(database), subtype=subname) rownames(datTraits)=datTraits[,1] head(datTraits) # Choose a set of soft-thresholding powers powers = c(c(1:10), seq(from = 12, to=20, by=2)) datExpr <- WGCNA_matrix # Call…

  • MA图绘制

    # 导入ggplot2包 library(ggplot2) # 设置好工作目录(到数据所在目录) # 读取输入数据 “R0-vs-R3.isoforms.filter.tsv” data = read.table(“R0-vs-R3.isoforms.filter.tsv”,header=T,row.names=1) # 计算M值和A值,并将M作为y轴,A作为x轴 aes = aes(x=(log2(R0_fpkm)+log2(R3_fpkm))/2,y=log2(R0_fpkm)-log2(R3_fpkm)) # 绘制MAplot ggplot(data=data,aes) + geom_point(aes(color=significant)) # 添加辅助线 ggplot(data=data,aes) + geom_hline(yintersect=0,linetype=4,color=”blue”) + geom_point(aes(color=significant)) # 改变点的大小 maplot = ggplot(data=data,aes) + geom_hline(yintercept=0,linetype=4,color=”blue”) + geom_point(aes(color=significant),size=1) # 设置自定义染色 maplot + scale_color_manual(values=c(“green”,”black”,”red”)) # 设置标题 maplot + scale_color_manual(values=c(“green”,”black”,”red”)) + labs(title=”MAplot of R0-vs-R3″,x=”A”,y=”M”) # 保存MAplot ggsave(“R0-vs-R3.MAplot.png”,width=8,height=6)

  • DESeq2使用流程

    library(“DESeq2”) database

  • PCA图绘制

    输入表达矩阵数据文件 library(ggplot2) library(gmodels) inname = “COS_deseq_counts_normalized.txt” outname = “COS_PCA.png” group

  • 绘制盒形图

    输入表达矩阵数据文件 library(ggplot2) library(reshape2) exprSet <- read.table(file = “COS_deseq_counts_normalized.txt”, sep = “\t”, header = T) group_list <- factor(c(rep(“COS0”,3),rep(“COS1”,3),rep(“COS3”,3),rep(“COS6”,3),rep(“COS12”,3)), levels = c(“COS0”, “COS1″,”COS3″,”COS6″,”COS12″)) exprSet_L=melt(exprSet) colnames(exprSet_L)=c(‘Gene_ID’,’sample’,’value’) exprSet_L$group=rep(group_list,each=nrow(exprSet)) #以sample作为x轴,group作为填充颜色 ggplot(exprSet,aes(x=sample,y=value,fill=group)) + geom_boxplot() #以group作为x轴,sample作为填充颜色 #ggplot(exprSet,aes(x=group,y=value,fill=sample)) + geom_boxplot() # 修改sample名称 exprSet_L$sample=paste(exprSet_L$group,exprSet_L$sample,sep=”-“) ggplot(exprSet_L,aes(x=sample,y=value,fill=group)) + geom_boxplot() # 离群点的大小 cvbox=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot(outlier.size=0) #outlier.size=-1 ,点将彻底消失 # 修改x轴字符方向 cvbox+theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1)) # 修改图例位置和主题 cvbox+theme_bw()+theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),legend.position=”top”) # geom_jitter():为图形设置抖动,防止相同点重叠 cvbox+theme_bw()+theme(axis.text.x=element_text(angle=90,vjust=0.5,hjust=1),legend.position=”top”) + geom_jitter()

  • 叔本华:作为意志和表象的世界(一)

    第一篇  世界作为表象初论  服从充分根据律的表象 经验和科学的客体 跳出童年时代吧,朋友,觉醒呵! ——J·J·卢梭 “世界是我的表象”:这是一个真理,是对于任何一个生活着和认识着的生物都有效的真理;不过只有人能够将它纳入反省的,抽象的意识罢了。并且,要是人真的这样做了,那么,在他那儿就出现了哲学的思考。于是,他就会清楚而确切地明白,他不认识什么太阳,什么地球,而永远只是眼睛,是眼睛看见太阳;永远只是手,是手感触着地球;就会明白围绕着他的这世界只是作为表象而存在着的,也就是说这世界的存在完全只是就它对一个其他事物的,一个进行“表象者”的关系来说的。这个进行“表象者”就是人自己。如果有一真理可以先验地说将出来,那就是这一真理了,因为这真理就是一切可能的、可想得到的经验所同具的那一形式的陈述。它比一切,比时间、空间、因果性等更为普遍,因为所有这些都要以这一真理为前提。我们既已把这些形式都认作根据律的一些特殊构成形态,如果其中每一形式只是对一特殊类型的表象有效,那么,与此相反,客体和主体的分立则是所有那些类型的共同形式。客体主体分立是这样一个形式:任何一个表象,不论是哪一种,抽象的或直观的,纯粹的或经验的,都只有在这一共同形式下,根本才有可能,才可想象。因此,再没有一个比这更确切,更不依赖其他真理,更不需要一个证明的真理了;即是说:对于“认识”而存在着的一切,也就是全世界,都只是同主体相关联着的客体,直观者的直观;一句话,都只是表象。当然,这里所说的对于现在,也对于任何过去,任何将来,对于最远的和近的都有效;因为这里所说的对于时间和空间本身就有效;而又只有在时间、空间中,所有这些[过去、现在、未来、远和近]才能区别出来。一切一切,凡已属于和能属于这世界的一切,都无可避免地带有以主体为条件[的性质],并且也仅仅只是为主体而存在。世界即是表象。 这个真理决不新颖。它已包含在笛卡儿所从出发的怀疑论观点中。不过贝克莱是断然把它说出来的第一人;尽管他那哲学的其余部分站不住脚,在这一点上,他却为哲学作出了不朽的贡献。康德首先一个缺点就是对这一命题的忽略,这在本书附录中将有详尽的交代。与此相反,吠檀多哲学被认为是毗耶舍的作品,这里所谈的基本原理在那里就已作为根本命题出现了,因此印度智者们很早就认识这一真理了。威廉·琼斯在他最近《论亚洲哲学》(《亚洲研究》,第四卷第164页)一文中为此作了证,他说:“吠檀多学派的基本教义不在于否认物质的存在,不在否认它的坚实性、不可入性、广延的形状(否认这些,将意味着疯狂),而是在于纠正世俗对于物质的观念,在于主张物质没有独立于心的知觉以外的本质,主张存在和可知觉性是可以互相换用的术语。”这些话已充分地表出了经验的实在性和先验的观念性两者的共存。 在这第一篇里,我们只从上述的这一方面,即仅仅是作为表象的一面来考察这世界。至于这一考察,虽无损于其为真理,究竟是片面的,从而也是由于某种任意的抽象作用引出来的,它宣告了每一个人内心的矛盾,他带着这一矛盾去假定这世界只是他的表象,另一方面他又再也不能摆脱这一假定。不过这一考察的片面性就会从下一篇得到补充,由另一真理得到补充。这一真理,可不如我们这里所从出发的那一个,是那么直接明确的,而是只有通过更深入的探讨,更艰难的抽象和“别异综同”的功夫才能达到的。它必然是很严肃的,对于每一个人纵不是可怕的,也必然是要加以郑重考虑的。这另一真理就是每人,他自己也能说并且必须说的:“世界是我的意志。” 在作这个补充之前,也就是在这第一篇里,我们必须坚定不移地考察世界的这一面,即我们所从出发的一面,“可知性”的一面:因此,也必须毫无抵触心情地将当前现成的客体,甚至自己的身体(我们就要进一步谈到这点)都仅仅作为表象看,并且也仅仅称之为表象。我们希望往后每一个人都会确切明白我们在这样做的时候,只仅仅是撇开了意志;而意志就是单独构成世界另外那一面的东西;因为这世界的一面自始至终是表象,正如另一面自始至终是意志。至于说有一种实在,并不是这两者中的任何一个方面,而是一个自在的客体(康德的“自在之物”可惜也不知不觉的蜕化为这样的客体),那是梦呓中的怪物;而承认这种怪物就会是哲学里引人误入迷途的鬼火。 转载:叔本华哲学智慧

  • GSEA命令行

    前边已经说到MSigDB数据库的entrez ID和symble ID都是人源的,其他物种要想做GSEA的话,就必须要ID转换,做成属于自己研究物种的数据库,这一步是关键。想要根据简单的Nr注释或其他注释来转换是行不通的,比如我研究的物种为硬骨鱼类,Nr注释结果分值最高的前20个一般也都是鱼类,就算注释到了人源的基因,也可能存在版本不对的问题。因此我们要知道MSigDB数据库基因ID的版本,然后下载对应的全基因组蛋白序列,然后比对卡阈值,做成自己的数据库。 GSEA的运行可以界面操作和命令行操作,界面操作教程很多就不多做赘述,它的不足之处除了操作繁琐之外,里面几个较大的库需要的内存较高,像我的4G渣渣电脑,好几个都运行失败,不得不转到服务器上运行;而且界面也不能批量运行。 java -cp gsea-3.0.jar -Xmx10000m xtools.gsea.Gsea -res ma/normalized_counts_all.gct -cls ma/samples.cls#COS1_versus_COS0 -gmx ma/ma.h.all.v6.1.symbols.gmt -collapse false -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_type no_balance -scoring_scheme weighted -rpt_label COS1vsCOS0_h -metric Ratio_of_Classes -sort real -order descending -create_gcts false -create_svgs false -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -rnd_seed timestamp -save_rnd_lists false…