面试了57位程序员，我发现混得好的人都有这个特质！

其实就是有解决问题的能力和意愿！与编程学习的宗旨不谋而合

没有严格意义上的真正编程能力，只有解决问题的能力。你解决的任何问题都依赖于别人解决过的子问题，所以不存在绝对的原创还是复用。

真正的编程能力其实并不是对语法细节的理解，也不在于手写或者复制粘贴，更不在于对什么操作系统的使用，或者常用库的api的记忆。而是找出解决方法的能力，把现实问题转换为代码逻辑的能力。这个是最重要的。语法很好学，只要看一看，再不行网上搜一搜都有，但是解决问题的能力，在网上搜不到，找不来，谁也帮不了。只能在长期的分析问题解决问题的过程中得到。

编程能力是一种解决问题的能力。如果问题没能被很好地解决，知道再多也没用。编程能力是一种运用机器解决问题的能力。首先是要判断问题在什么程度上可被机器解决，比如理论计算机科学会告诉我们什么可做、什么理论上不可做、什么理论上可做实践上不可做。然后是让机器更好地理解问题，比如计算机都是（图灵－冯诺依曼模型）等价，但不同的问题可能会适用不同的编程语言。再后是让机器能更高效率地解决问题，比如同样的问题可能会有效率差别巨大的算法。

编程能力是一种预见未来的能力。目前的方案有哪些假设和局限性，在何种情形下会遇到问题甚至崩溃。在未来出现问题时问题是否需要重新定义，系统是否需要重新设计，代码是否需要重构或优化等等都需要未雨绸缪。

其实不管是复制黏贴也好，自己手写也好，关键的是解决问题。编程最终还是个生产工具，目的是解决问题，不能解决问题的，一切都是空中楼阁，毫无价值。

下面是一个解决问题的示例
对一个seurat的对象 里面的全部单细胞进行gsva分析会面临内存问题

如果跑gsva面临内存问题

代码如下：

sce=pbmc3k
mat=sce@assays$RNA@counts
mat[1:4,1:4]
X=as.matrix(mat)
dim(X)

gmtfile ='../MSigDB/symbols/h.all.v7.2.symbols.gmt'
library(GSEABase)
library(GSVA)
geneset <- getGmt( gmtfile )  
geneset
# method=c("gsva", "ssgsea", "zscore", "plage"),
# kcdf=c("Gaussian", "Poisson", "none"),
# mx.diff=FALSE: ES is calculated as the maximum distance of the random walk from 0. 
# mx.diff=TRUE (default): ES is calculated as the magnitude difference between the largest positive and negative random walk deviations.
es.max <- gsva(X, geneset, 
   mx.diff=FALSE, verbose=FALSE, 
   parallel.sz= 8 )

其中 h.all.v7.2.symbols.gmt 文件需要自行下载哦，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

  • H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
  • C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
  • C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
  • C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
  • C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
  • C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
  • C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO  发表芯片数据
  • C7: immunologic signatures: 免疫相关基因集合。

因为这个案例，仅仅是3千个单细胞，所以很难遇到内存不足的报错，其实真实情况下 往往是5万个单细胞以上，大家如果有自己的数据，就可以测试了。

其实解决内存不足问题超级简单

只需要把前面的单细胞的表达矩阵，根据细胞类型拆分开来，分开独立跑gsva，然后再合并，代码如下：

table(sce@meta.data$orig.ident)
es.max.list <- lapply(unique(sce@meta.data$orig.ident),function(i){
  Y=X[,sce@meta.data$orig.ident == i ]
  es.max <- gsva(Y, geneset, 
  mx.diff=FALSE, verbose=FALSE, 
  parallel.sz=8)
  return(es.max)
}) 
es.max=do.call(cbind,es.max.list)
es.max[1:4,1:4]

当然了，这个时候表达矩阵的顺序其实是发生了变化哦，后续如果对这个gsva矩阵绘制热图的时候，需要调整回来哦。

es.max[1:4,1:4]
X[1:4,1:4]
identical(colnames(es.max),colnames(X))
pos=match(colnames(X),colnames(es.max))
es.max=es.max[,pos]
identical(colnames(es.max),colnames(X))

其实呢，也不一定要根据细胞类型来拆分表达量矩阵分开跑gsva，只需要拆分即可。毕竟啊， 按照细胞类型来拆分表达量矩阵，通常呢，并不是等价的。

 as.data.frame(table(pbmc3k$seurat_annotations))
 Var1 Freq
1  Naive CD4 T  697
2 Memory CD4 T  483
3   CD14+ Mono  480
4   B  344
5  CD8 T  271
6 FCGR3A+ Mono  162
7  NK  155
8  DC   32
9  Platelet   14

最后的可视化如下

ac=data.frame(celltype=pbmc3k$seurat_annotations)
rownames(ac)=rownames(sce@meta.data)
pheatmap::pheatmap(es.max,
 annotation_col = ac,
 show_colnames = F)

挺漂亮的：

生信小技巧
注意看 https://www.bilibili.com/video/BV1hW411c7Zi (复制粘贴到浏览器打开)

理论上我推荐大家一定学会R语言的，本次视频我虽然讲解的是perl，但实际上是在讲解编程思维，请仔细听P6:

  • 视频里面讲解为什么学习编程

  • 一个例子

  • 编程思维的获得

    • perl -alne '{if($F[2] eq "transcript" ){/transcript_id\s\"(.*?)\"/;$transcript_id=$1; /FPKM\s\"(.*?)\"/;$fpkm=$1; /cov\s\"(.*?)\"/;$cov=$1; print "$transcript_id\t$fpkm\t$cov" } }' A172D2_1.gtf|head

    • perl -alne '{if($F[2] eq "transcript" ){/transcript_id\s\"(.*?)\"/;$transcript_id=$1; /FPKM\s\"(.*?)\"/;$fpkm=$1; /cov\s\"(.*?)\"/;$cov=$1;  } else{$c{$transcript_id}++;$l{$transcript_id}+=(abs($F[4]-$F[3])+1)} } END{print "$_\t$c{$_}\t$l{$_}" foreach keys %c}' A172D2_1.gtf|grep ENST00000446136.1

    • perl -alne '{if($F[2] eq "transcript" ){/transcript_id\s\"(.*?)\"/;$transcript_id=$1; /FPKM\s\"(.*?)\"/;$fpkm=$1; /cov\s\"(.*?)\"/;$cov=$1; $f{$transcript_id}=$fpkm;$c{$transcript_id}=$cov } else{$c{$transcript_id}++;$l{$transcript_id}+=(abs($F[4]-$F[3])+1)} } END{print "$_\t$c{$_}\t$l{$_}\t$f{$_}\t$c{$_}" foreach keys %c}' A172D2_1.gtf|grep ENST00000446136.1

    • perl -alne '{if($F[2] eq "transcript" ){/transcript_id\s\"(.*?)\"/;$transcript_id=$1; /FPKM\s\"(.*?)\"/;$fpkm=$1; /cov\s\"(.*?)\"/;$cov=$1; $f{$transcript_id}=$fpkm;$c{$transcript_id}=$cov } else{$counts{$transcript_id}++;$l{$transcript_id}+=(abs($F[4]-$F[3])+1)} } END{foreach (keys %c) {next if $f{$_} <1;next if $c{$_} <3; next if $counts{$_} <2;next if $l{$_} <200;print "$_\t$counts{$_}\t$l{$_}\t$f{$_}\t$c{$_}" }}' A172D2_1.gtf  >need.transcripts

    • 完美解决问题

    • 读取文件

    • 默认变量，关键词

    • 按列拆分

    • 判断语句

    • 循环语句

    • ①FPKM>=1;

    • ②cov >=3;

    • ③exon>=2；

    • ④exon总长度>=200

    • 解析需求

    • 熟悉语法

    • 不停的调试

  • 学什么语言

    • awk 建议了解基础

    • perl 建议不学

    • python 强烈推荐 （小学生都会）

    • java 等其它可以了解（方便转行）

  • 学习资源：纸质版书籍 https://mp.weixin.qq.com/s/Gks3ogPPpsagocxjeBp8dQ

  • 练习为王：生信编程实战；https://ke.qq.com/course/285055

  • 生信编程200题：http://www.biotrainee.com/forum-90-1.html

赞 (0)

相关推荐