十 30

R里面list对象和AnnDbBimap对象区别

Posted on by

AnnDbBimap对象是R里面的bioconductor系列包的基础对象，在探针数据里面会包装成ProbeAnnDbMap，跟go通路相关又是GOTermsAnnDbBimap对象。

但是它都是AnnDbBimap对象衍生过来的

主要存在于芯片系列的包和org系列的包，其实AnnDbBimap对象就是list对象的包装，比如下面这些例子：

ls("package:hgu133plus2.db")

 [1] "hgu133plus2"              "hgu133plus2.db"           "hgu133plus2_dbconn"
 [4] "hgu133plus2_dbfile"       "hgu133plus2_dbInfo"       "hgu133plus2_dbschema"
 [7] "hgu133plus2ACCNUM"        "hgu133plus2ALIAS2PROBE"   "hgu133plus2CHR"
[10] "hgu133plus2CHRLENGTHS"    "hgu133plus2CHRLOC"        "hgu133plus2CHRLOCEND"
[13] "hgu133plus2ENSEMBL"       "hgu133plus2ENSEMBL2PROBE" "hgu133plus2ENTREZID"
[16] "hgu133plus2ENZYME"        "hgu133plus2ENZYME2PROBE"  "hgu133plus2GENENAME"
[19] "hgu133plus2GO"            "hgu133plus2GO2ALLPROBES"  "hgu133plus2GO2PROBE"
[22] "hgu133plus2MAP"           "hgu133plus2MAPCOUNTS"     "hgu133plus2OMIM"
[25] "hgu133plus2ORGANISM"      "hgu133plus2ORGPKG"        "hgu133plus2PATH"
[28] "hgu133plus2PATH2PROBE"    "hgu133plus2PFAM"          "hgu133plus2PMID"
[31] "hgu133plus2PMID2PROBE"    "hgu133plus2PROSITE"       "hgu133plus2REFSEQ"
[34] "hgu133plus2SYMBOL"        "hgu133plus2UNIGENE"       "hgu133plus2UNIPROT"

那么我们可以随便挑选包中的一个数据分析一下

x <- hgu133plus2SYMBOL

xlist=as.list(x)

我们查看X对象，可知，它是object of class "ProbeAnnDbMap"，而这个对象，就是常见的list对象，被包装了一下， 只有我们明白了它和list对象到底有什么区别，才算是真正搞懂了这一系列包

但是这个ProbeAnnDbMap对象，在GO等包里面会被包装成更复杂的对象-GOTermsAnnDbBimap，但是对他们的理解都大同小异。

我们先回顾一下在R语言里面的list的基础知识：

R 的 列表（list）是一个以对象的有序集合构成的对象。 列表中包含的对象又称为它的分量（components）。

分量可以是不同的模式或者类型， 如一个列表可以包括数值向量，逻辑向量，矩阵，复向量， 字符数组，函数等等

如果你会perl的话，可以把它理解成hash。

分量常常会被编号的（numbered），并且可以利用它来 访问分量

列表的分量可以被命名，这种情况下 可以通过名字访问。

特别要注意一下 Lst[[1]] 和 Lst[1] 的差别。 [[...]] 是用来选择单个 元素的操作符，而 [...] 是一个更为一般 的下标操作符。

因此前者得到的是列表 Lst 中的第一个对象, 并且含有分量名字的命名列表（named list）中分量的名字会被排除在外的。

后者得到的则是列表 Lst 中仅仅由第一个元素 构成的子列表。如果是命名列表， 分量名字会传给子列表的。

那么接下来我们就看看x和xlist的区别。

它们里面的数据都是affymetrix公司出品的人类的hgu133plus2芯片的探针ID与基因symbol的对应关系

如果我想拿到所有的探针ID，那么对于AnnDbBimap对象，需要用mappedkeys(x)，对于普通的list对象，需要用names(xlist).

对于普通的list对象，如果我想看前几个元素，直接head就可以了，但是对于AnnDbBimap对象，数据被封装好了，需要先as.list，然后才能head

mapped_probes <- mappedkeys(x)

PID2=head(mapped_probes)

那么，如果我们想根据以下探针ID来查看它们在这些数据里面被对应着哪些基因symbol 呢？

> PID2 #这是一串探针ID，后面的操作都需要用的

[1] "1053_at"   "117_at"    "121_at"    "1255_g_at" "1316_at"   "1320_at"

如果是对于AnnDbBimap对象，我们可以用mget函数来操作，取多个探针的对应基因symbol

accnum <- mget(PID2, env=hgu133plus2ACCNUM);

gname <- mget(PID2, env=hgu133plus2GENENAME)

gsymbol <- mget(PID2, env=hgu133plus2SYMBOL)

mget函数返回的就是普通的list函数了，可以直接查看了。

如果是对于普通的list对象，我们想取多个探针的对应基因symbol也是非常简单的，xlist[PID2]即可。

那么我们不禁有问了，既然它们两个功能完全一样，何苦把它包装成一个对象了，我直接操作list对象不就好了，学那么多规矩干嘛？

所以，重点就来了

>  length(mapped_probes)

[1] 42125

> length(names(xlist))

[1] 54675

看懂了吗？

 

但是，事实上用处也不大，你觉得下面这两个有区别吗？

SYMBOL <- AnnotationDbi::as.list(hgu133plus2SYMBOL)

SYMBOL <- SYMBOL[!is.na(SYMBOL)];

 

x <- hgu133plus2SYMBOL

mapped_probes <- mappedkeys(x)

SYMBOL <- AnnotationDbi::as.list(x[mapped_probes])

 

PS，在R里面创建一个list对象是非常简单的

The setNames() built-in function makes it easy to create a hash from given key and value lists.(Thanks to Nick K for the better suggestion.)

Usage: hh <- setNames(as.list(values), keys)

Example:

players <- c("bob", "tom", "tim", "tony", "tiny", "hubert", "herbert")
rankings <- c(0.2027, 0.2187, 0.0378, 0.3334, 0.0161, 0.0555, 0.1357)
league <- setNames(as.list(rankings), players)

十 29

gene的各种ID转换终结者-bioconductor系列包

Posted on by

经常会有人问这样的问题I have list of 10,000 Entrez IDs and i want to convert the multiple Entrez IDs into the respective gene names. Could someone suggest me the way to do this?等等类似的基因转换，能做的基因转换的方法非常多，以前不懂编程的时候，都是用各种网站，而最常用的就是ensembl的biomart了，它支持的ID非常多，高达几百种，好多ID我到现在都不知道是什么意思。

现在学会编程了，我比较喜欢的是R的一些包，是bioconductor系列，一般来说，其中有biomart，org.Hs.eg.db，annotate，等等。关于biomart我就不再讲了，我前面的博客至少有七八篇都提到了它。本次我们讲讲简单的， 我就以把gene entrez ID转换为gene symbol 为例子把。

当然，首先要安装这些包，并且加载。

if("org.Hs.eg.db" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("org.Hs.eg.db")}
suppressMessages(library(org.Hs.eg.db))  我比较喜欢这样加载包

library(annotate) #一般都是这样加载包

如果是用org.Hs.eg.db包，首先你只需要读取你的待转换ID文件，构造成一个向量，tmp，然后只需要symbols <- org.Hs.egSYMBOL[as.character(tmp)]就可以得到结果了，返回的symbols是一个对象，需要用toTable这个函数变成数据框。但是这样转换容易出一些问题，比如如果你的输入数据tmp，里面含有一些无法转换的gene entrez ID，就会报错。

而且它支持的ID转换很有限，具体看看它的说明书即可：https://www.bioconductor.org/packages/release/data/annotation/manuals/org.Hs.eg.db/man/org.Hs.eg.db.pdf

org.Hs.eg.db
org.Hs.eg_dbconn
org.Hs.egACCNUM
org.Hs.egALIAS2EG
org.Hs.egCHR
org.Hs.egCHRLENGTHS
org.Hs.egCHRLOC
org.Hs.egENSEMBL
org.Hs.egENSEMBLPROT
org.Hs.egENSEMBLTRANS
org.Hs.egENZYME
org.Hs.egGENENAME
org.Hs.egGO
org.Hs.egMAP
org.Hs.egMAPCOUNTS
org.Hs.egOMIM
org.Hs.egORGANISM
org.Hs.egPATH
org.Hs.egPMID
org.Hs.egREFSEQ
org.Hs.egSYMBOL
org.Hs.egUCSCKG
org.Hs.egUNIGENE
org.Hs.egUNIPROT

如果是用annotate包，首先你还是需要读取你的待转换ID文件，构造成一个向量，tmp，然后用getSYMBOL(as.character(tmp), data='org.Hs.eg')这样直接就返回的还是以向量，只是在原来向量的基础上面加上了names属性。说明书：http://www.bioconductor.org/packages/3.3/bioc/manuals/annotate/man/annotate.pdf

然后你可以把转换好的向量写出去，如下：

1 A1BG
2 A2M
3 A2MP1
9 NAT1
10 NAT2
12 SERPINA3
13 AADAC
14 AAMP
15 AANAT
16 AARS

PS：如果是芯片数据，需要把探针的ID转换成gene，那么一般还需要加载特定芯片的数据包才行：

platformDB <- paste(eset.mas5@annotation, ".db", sep="") #这里需要确定你用的是什么芯片
cat("the annotation is ",platformDB,"\n")
if(platformDB %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");tmp=try(biocLite(platformDB))}
library(platformDB, character.only=TRUE)
probeset <- featureNames(eset.mas5)
rowMeans <- rowMeans(exprSet)

library(annotate) # lookUp函数是属于annotate这个包的
EGID <- as.numeric(lookUp(probeset, platformDB, "ENTREZID"))

十 29

使用GEOmetadb包来获取对应GEO数据的实验信息

Posted on by

理论上我前面提到的GEOquery包就可以根据一个GSE索引号来获取NCBI提供的所有关于这个GSE索引号的数据了，包括metadata，表达矩阵，soft文件，还有raw data

但是很多时候，那个metadata并不是很整齐，而且一个个下载太麻烦了，所以就需要用R的bioconductor的另一个神奇的包了GEOmetadb

它的示例：http://bioconductor.org/packages/devel/bioc/vignettes/GEOmetadb/inst/doc/GEOmetadb.R

里面还是很多数据库基础知识的
代码托管在github，它的示例代码是这样连接数据库的：
library(GEOmetadb)
if(!file.exists('GEOmetadb.sqlite')) getSQLiteFile()
file.info('GEOmetadb.sqlite')
con <- dbConnect(SQLite(),'GEOmetadb.sqlite')
dbDisconnect(con)
但是一般不会成功，因为这个包把它的GEOmetadb.sqlite文件放在了国外网盘共享，在国内很难访问，推荐大家想办法下载到本地
tmp2
用这个代码就会成功了，需要自己下载GEOmetadb.sqlite文件然后放在指定目录：/path/GEOmetadb.sqlite 需要自己修改
我们的diabetes.GEO.list文件内容如下：
GSE1009
GSE10785
GSE1133
GSE11975
GSE121
GSE12409
那么会产生的表格文件如下：共有32列数据信息，算是蛮全面的了
tmp
十 28

模拟Y染色体测序判断，并比对到X染色体上面，看同源性

Posted on by

首先下载两条染色体序列

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz;

wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz;

152M Mar 21  2009 chrX.fa

58M Mar 21  2009 chrY.fa

然后把X染色体构建bwa的索引

bwa index chrX.fa

[bwa_index] Pack FASTA... 1.97 sec

[bwa_index] Construct BWT for the packed sequence...

[BWTIncCreate] textLength=310541120, availableWord=33850812

[BWTIncConstructFromPacked] 10 iterations done. 55838672 characters processed.

[BWTIncConstructFromPacked] 20 iterations done. 103157920 characters processed.

[BWTIncConstructFromPacked] 30 iterations done. 145211344 characters processed.

[BWTIncConstructFromPacked] 40 iterations done. 182584528 characters processed.

[BWTIncConstructFromPacked] 50 iterations done. 215797872 characters processed.

[BWTIncConstructFromPacked] 60 iterations done. 245313968 characters processed.

[BWTIncConstructFromPacked] 70 iterations done. 271543920 characters processed.

[BWTIncConstructFromPacked] 80 iterations done. 294853104 characters processed.

[bwt_gen] Finished constructing BWT in 88 iterations.

[bwa_index] 98.58 seconds elapse.

[bwa_index] Update BWT... 0.96 sec

[bwa_index] Pack forward-only FASTA... 0.91 sec

[bwa_index] Construct SA from BWT and Occ... 33.18 sec

[main] Version: 0.7.8-r455

[main] CMD: /lrlhps/apps/bioinfo/bwa/bwa-0.7.8/bwa index chrX.fa

[main] Real time: 141.623 sec; CPU: 135.605 sec

由于X染色体也就152M，所以很快，两分钟解决战斗！

然后模拟Y染色体的测序判断（PE100，insert400）

209M Oct 28 10:19 read1.fa

209M Oct 28 10:19 read2.fa

模拟的程序很简单

tmp

 

while(<>){
chomp;
$chrY.=uc $_;
}
$j=0;
open FH_L,">read1.fa";
open FH_R,">read2.fa";
foreach (1..4){
for ($i=600;$i<(length($chrY)-600);$i = $i+50+int(rand(10))){
$up = substr($chrY,$i,100);
$down=substr($chrY,$i+400,100);
next unless $up=~/[ATCG]/;
next unless $down=~/[ATCG]/;
$down=reverse $down;
$down=~tr/ATCG/TAGC/;
$j++;
print FH_L ">read_$j/1\n";
print FH_L "$up\n";
print FH_R ">read_$j/2\n";
print FH_R "$down\n";
}
}
close FH_L;
close FH_R;

然后用bwa mem 来比对

bwa mem -t 12 -M chrX.fa read*.fa >read.sam

用了12个线层，所以也非常快

[main] Version: 0.7.8-r455

[main] CMD: /apps/bioinfo/bwa/bwa-0.7.8/bwa mem -t 12 -M chrX.fa read1.fa read2.fa

[main] Real time: 136.641 sec; CPU: 1525.360 sec

643M Oct 28 10:24 read.sam

然后统计比对结果

samtools view -bS read.sam >read.bam

158M Oct 28 10:26 read.bam

samtools flagstat read.bam

3801483 + 0 in total (QC-passed reads + QC-failed reads)

0 + 0 duplicates

2153410 + 0 mapped (56.65%:-nan%)

3801483 + 0 paired in sequencing

1900666 + 0 read1

1900817 + 0 read2

645876 + 0 properly paired (16.99%:-nan%)

1780930 + 0 with itself and mate mapped

372480 + 0 singletons (9.80%:-nan%)

0 + 0 with mate mapped to a different chr

0 + 0 with mate mapped to a different chr (mapQ>=5)

我自己看sam文件也发现真的同源性好高呀，总共就模拟了380万reads，就有120万是百分百比对上了。

所以对女性个体来说，测序判断比对到Y染色体是再正常不过的了。如果要判断性别，必须要找那些X,Y差异性区段

十 23

生物信息小白如何自学编程

Posted on by

这本来是我在知乎上面看到的问题，所以就抽空回答了一下：http://www.zhihu.com/question/36701137/answer/68928111

首先，你懂得想去看源码，这是一个很好的兆头，一些非常正规的源码的确是编程进阶的的捷径，毕竟我们大部分人都不可能得到别人的手把手指导，所以只能靠自己的悟性了。

我就以我自己的经历来回答这个问题吧，我作为一个纯生物出身的小白，现在编程技术应该还算可以了！

首先，不管是哪个语言，perl,python,R,matlab都好，它们都有一堆的基础书籍，你必须以囫囵吞枣的心态看完一两本书（书没有好坏，别要我给你推荐书名），必须看完，了解编程基础。

接下来的步骤最重要，就是实践，不停的实践，在实践中运用编程技术，这样是学的最快的，不然你看再多的书也只是一个概念。

我这里重点推荐一个工具集，它实现了很多生物信息学需要的常用操作，网址是：Bioinformatics Tools
包含以下64中工具，而且网页也很清楚的描述了它们的功能，其实非常简单，但是这样写程序非常有效。
"Combines multiple FASTA entries into a single sequence."
"Returns the entire sequence contained in an EMBL file in FASTA format."
"Parses the feature table of an EMBL file and returns the feature sequences."
"Parses the feature table of an EMBL file and returns the protein translations."
"Removes non-DNA characters from text."
"Removes non-protein characters from text."
"Returns the entire sequence contained in a GenBank file in FASTA format."
"Parses the feature table of a GenBank file and returns the feature sequences."
"Parses the feature table of a GenBank file and returns the protein translations."
"Converts single letter amino acid codes to three letter codes."
"Reads a list of positions and ranges and returns those parts of a DNA sequence."
"Reads a list of positions and ranges and returns those parts of a protein sequence."
"Determines the reverse-complement, reverse, or complement of the sequence you enter."
"Separates bases according to codon position."
"Converts a FASTA sequence into multiple sequences."
"Converts three letter amino acid codes to one letter codes."
"Returns DNA sequence segments specified by a position and window size."
"Returns protein sequence segments specified by a position and window size."
"Plots codon frequency (according to the codon table you enter) for each codon in a DNA sequence."
"Returns a standard codon usage table."
"Returns a list of potential CpG islands."
"Calculates the molecular weight of DNA sequences."
"Returns positions of the patterns you enter."
"Returns basic sequence statistics."
"Returns sequences that are identical or similar to a query sequence."
"Returns sequences that are identical or similar to a query sequence."
"Accepts aligned sequences in FASTA format and calculates the identity and similarity of each sequence pair."
"Can be used to predict a DNA sequence in another species using a protein sequence alignment."
"Finds DNA sequences that can easily be converted to a restriction site."
"Determines the positions of open reading frames."
"Returns the optimal global alignment for two coding DNA sequences."
"Returns the optimal global alignment for two DNA sequences."
"Returns the optimal global alignment for two protein sequences."
"Returns a report describing PCR primer properties"
"Generates PCR products from a template and two primer sequences."
"Returns the grand average of hydropathy value of protein sequences."
"Returns the predicted isoelectric point of protein sequences."
"Calculates the molecular weight of protein sequences."
"Returns positions of the patterns you enter."
"Returns basic sequence statistics."
"Converts the sequence you enter into restriction fragments."
"Returns the number and positions of restriction sites."
"Can be used to convert protein into DNA."
"Returns the translation in the reading frame you specify."
"Colors a sequence alignment based on sequence conservation."
"Colors a protein alignment based on biochemical properties of residues."
"Numbers and groups DNA according to your specifications."
"Numbers and groups amino acids according to your specifications."
"Shows PCR primer annealing sites, translations, and restriction sites."
"Shows restriction sites and protein translations."
"Shows protein translations."
"Introduces random mutations into DNA sequences."
"Introduces random mutations into protein sequences."
"Generates a random coding sequence of the length you specify."
"Generates a random DNA sequence of the length you specify."
"Replaces regions of the DNA sequences you enter with random bases."
"Generates a random protein sequence of the length you specify."
"Replaces regions of the protein sequences you enter with random residues."
"Samples bases from a DNA sequence with replacement."
"Samples residues from a protein sequence with replacement."
"Randomly shuffles the DNA sequences you enter."
"Randomly shuffles the protein sequences you enter."
"IUPAC codes for DNA and protein."
"The genetic codes used in the Sequence Manipulation Suite."
当你实现完了这些需求，你不仅仅学会了编程，而且是学会了编程该如何应用在生物信息学里面！
用perl,python,R,matlab中的任何一种都可以实现，它们没有任何区别的，别纠结语言的问题。
不推荐初学者看源代码，因为源代码太正规了，定义变量就几十行代码了，再定义函数又是几百行代码，而真正学生物信息学的压根写代码都不超过五十行的，比如我上面提到那64个生物数据处理需求，一般就七八行代码就可以（在perl里面）
不信你可以看看这个github里面托管的代码：trinityrnaseq/util/misc at master · trinityrnaseq/trinityrnaseq · GitHub
里面有很多perl代码，都是实现各种数据转换的，写的非常正规，甚至能把一行代码就能解决的问题写成几百甚至上千行，除非你想把自己的代码拿去发文章或者出售，否则正常的生物信息学研究根本用不着！
当然，回到你最初的问题，哪里能找到源码呢？
首先，你可以去图书馆看一堆书籍，它们都会有光盘，下载既有视频又有源码，或者书上一般会说源码在哪里下载，比如这个pleac/include/perl at master · pleac/pleac · GitHub
然后，你可以找一大堆的生物信息学软件，它们一般都托管在github上面，这个链接里面有三百多个生物信息学转录组领域的软件：List of RNA-Seq bioinformatics tools
这个链接有几百个生物信息学里面做alignment的软件：
甚至连常见的生物信息学数据库也有自己的源码包：例如NCBI，ensembl，UCSC
下面就是ENSEMBL数据库的：NGS数据比对工具持续收集
（记住，这些软件都是人家发表文章的，非常难，你一辈子能搞定一个就很了不起了，比如我，就搞了一下bowtie，也是一知半解的）
分享了所有的代码，实在是太方便了：Ensembl Project · GitHub
可以跟着这些代码学习编程：Ensembl/ensembl-pipeline · GitHub
它的官网的帮助文档也特别详细：Help & Documentation
你现在还缺资料吗？

十 22

R语言中的排序，集合运算，reshape，以及merge总结

Posted on by

一直以来我都是随便看了点R的编程教程，因为我学了一点点C，所以还算有基础，现在基本上简单看看教程就能懂一门语言了，区别只是熟练度而已。R用得比较多，所以还算擅长，但是很多快捷应用的地方，我总是寄希望于到时候再查资料，所以没能用心的记住，这次花了点时间好好整理了一下R里面关于数据操作的重点，我想，以后再碰到类似的数据处理要求，应该很快能解决了把。

首先看看排序：

在R中，和排序相关的函数主要有三个：sort()，rank()，order()。

sort(x)是对向量x进行排序，返回值排序后的数值向量。

rank()是求秩的函数，它的返回值是这个向量中对应元素的“排名”。

order()的返回值是对应“排名”的元素所在向量中的位置。

其中sort(x)等同于x[order(x)]

下面以一小段R代码来举例说明：

> x<-c(97,93,85,74,32,100,99,67)

> sort(x)

[1]  32  67  74  85  93  97  99 100

> order(x)

[1] 5 8 4 3 2 1 7 6

> rank(x)

[1] 6 5 4 3 1 8 7 2

> x[order(x)]

[1]  32  67  74  85  93  97  99 100

其中比较有用的order，它可以用来给数据框进行排序

dat[order(dat[,1]),] 以该数据框的第一列进行排序

dat[order(dat[,1],dat[,2]),] 以该数据框的第一列为主要次序，第二列为次要序列进行排序

 

然后我们看看集合运算：

在R里面除了简单的对两个向量求交集并集补集之外，比较重要的就是match和 %in% 了，需要重点讲讲。

#首先对集合A,B,C赋值

> A<-1:10

> B<-seq(5,15,2)

> C<-1:5

> #求A和B的并集

> union(A,B)

[1]  1  2  3  4  5  6  7  8  9 10 11 13 15

> #求A和B的交集

> intersect(A,B)

[1] 5 7 9

> #求A-B

> setdiff(A,B)

[1]  1  2  3  4  6  8 10

> #求B-A

> setdiff(B,A)

[1] 11 13 15

> #检验集合A,B是否相同

> setequal(A,B)

[1] FALSE

> #检验元素12是否属于集合C

> is.element(12,C)

[1] FALSE

> #检验集合A是否包含C

> all(C%in%A)

[1] TRUE

> all(C%in%B)

从上面可以看到%in%这个操作符只返回逻辑向量TRUE 或者FALSE，而且返回值应该与%in%这个操作符前面的向量程度相等。也就是说它相当于遍历了C里面的一个个元素，判断它们是否在B中出现过，然后返回是或者否即可。

而match(C,B)的结果就很不一样了，它的返回结果同样与前面的向量等长，但是它并非返回逻辑向量，而是遍历了C里面的一个个元素，判断它们是否在B中出现过，如果出现就返回在B中的索引号，如果没有出现，就返回NA。

>B

[1]  5  7  9 11 13 15

>C

[1] 1 2 3 4 5

>match(C,B)

[1] NA NA NA NA  1

>C%in%B

[1] FALSE FALSE FALSE FALSE  TRUE

 接下来我们看看reshape：

这是一个需要安装的包，起得就是R里面最经典的把长型数据变宽，和把宽数据拉长的作用。

其中melt函数是把很宽的数据拉长，它就是需要指定几列数据是保证不被融合的， 其余每一列数据都必须被融合到一列了，融合后的这一列数据每个元素旁边就用列名来标记该数据来自于哪一列。

# example of melt function
library(reshape)
mdata <- melt(mydata, id=c("id","time"))

融合后的数据如下：

id time variable value
1 1 x1 5
1 2 x1 3
2 1 x1 6
2 2 x1 2
1 1 x2 6
1 2 x2 5
2 1 x2 1
2 2 x2 4

可以看到variable列里面有两个不同的元素，说明是把旧数据中的两列给融合了，融合后的一个很长的列就是value

而cast函数的功能就是把刚才融合好的数据给还原。

# cast the melted data
# cast(data, formula, function)
subjmeans <- cast(mdata, id~variable, mean)
timemeans <- cast(mdata, time~variable, mean)

可以看到它们返回的结果是：

subjmeans

id x1 x2
1 4 5.5
2 4 2.5

timemeans

time x1 x2
1 5.5 3.5
2 2.5 4.5

可以看到cast函数比较复杂一点，formula公式右边的变量是需要拆开的variable，这一列变量有多少不重复元素，就新增多少列，左边的变量是行标记了，其余所有数据都会被计算一下再放在合适的位置。

在reshape2这个包里面还进化出来dcast和acast函数，功能大同小异。

 

最后我们来看看merge函数：

这个函数的功能非常强大，类似于SQL语句里面的join系列函数

测试数据如下，它们这两个表的连接是作者名

Capture1

 

如果要实现类似sql里面的inner join 功能，则用代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name")

如果要实现left join功能则用代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all.x=TRUE)

right join功能代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all.y=TRUE)

all join功能代码

m1 <- merge(authors, books, by.x = "surname", by.y = "name",all=TRUE)

关于单变量匹配的总结就是这些，但对于多变量匹配呢，例如下面两个表，需要对k1,k2两个变量都相等的情况下匹配

x <- data.frame(k1 = c(NA,NA,3,4,5), k2 = c(1,NA,NA,4,5), data = 1:5)

y <- data.frame(k1 = c(NA,2,NA,4,5), k2 = c(NA,NA,3,4,5), data = 1:5)

匹配代码如下merge(x, y, by = c("k1","k2"))  #inner join

另外一个多行匹配的例子如下：

Capture2

我们的测试数据如上，这两个表的连接在于作者名。下面这两个语句等价

merge(authors, books, by=1:2)

merge(authors, books, by.x=c("FirstName", "LastName"),

by.y=c("AuthorFirstName", "AuthorLastName"),

all.x=TRUE)

都可以把两张表关联起来。

当然，在我搜索资料的时候，发现了另外一个解决问题的方法：

A[with(A, paste(C1, C2, sep = "\r")) %in% with(B, paste(C1, C2, sep="\r")), ]

 

 

参考：http://blog.sina.com.cn/s/blog_6caea8bf0100spe9.html

http://blog.sina.com.cn/s/blog_6caea8bf010159dt.html

http://blog.csdn.net/u014801157/article/details/24372441

http://www.statmethods.net/management/reshape.html

https://docs.tibco.com/pub/enterprise-runtime-for-R/1.5.0_may_2013/TERR_1.5.0_LanguageRef/base/merge.html

http://r.789695.n4.nabble.com/Matching-multiple-columns-in-a-data-frame-td890832.html

http://bbs.pinggu.org/thread-3234639-1-1.html

http://www.360doc.com/content/14/0821/13/18951041_403561544.shtml

http://developer.51cto.com/art/201312/423612.htm

十 22

用R和mysql把基因对应的最大表达量探针挑出来

Posted on by

懂点芯片数据分析的都应该知道，芯片设计的时候，对一个基因设计了多个探针，这样设计是为了更好的捕获某些难以发现的基因，或者重点研究某些基因。
但是对我们的差异分析不方便，所以我们只分析哪些有对应了entrez ID的探针，而且对每个entrez ID，我们只需要挑选它表达量最高的那个探针。
所以就演化为一个编程问题：分组求最大值，多公共列合并！
如果是在R语言里面，那么首先这个table的表示形式如下
> esetDataTable[1:10,c(7,8)]
EGID rowMeans
1000_at 5595 1840.04259751826
1001_at 7075 799.075414422572
1002_f_at 1557 50.4884096416177
1003_s_at 643 142.372008051308
1004_at 643 211.65300963049
1005_at 1843 4281.29318032004
1006_at 4319 38.5784289213085
1007_s_at NA 1489.98158531843
1008_f_at 5610 4013.576753977
1009_at 3094 3070.50648167305
我们首先看看一个R语言函数处理方式吧，这个是比较容易想到的算法，但是用到了循环，非常的不经济，计算量很大。

Capture1
因为里面涉及到了双重循环。我进行了人工计时，这个程序耗费了一分钟十二秒，程序里面的计时器有点问题。

然后我再讲一个精简版的算法
dat=esetDataTable[,c(7,8)]
dat=as.data.frame(apply(dat,2,as.numeric))
dat$prob=rownames(esetDataTable)
首先可以得到需要提取的数据所在行的两个值
match_row=aggregate(dat,by=list(dat[,1]),max)
类似于这句SQL语句：SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID
现在要根据match_row表去原表esetDataTable里面提取我们的探针ID数据。

Capture2

这样就完美解决了问题，我们的合并后的数据的prob列，就是前面那个函数计算了一分多钟的返回的非冗余探针ID向量，relevantProbesets，但是这次只花了不到一秒钟。
tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans"))
setdiff(as.character(relevantProbesets),as.character(tmp_prob$prob))
length(union(as.character(relevantProbesets),as.character(tmp_prob$prob)))
setdiff(as.character(tmp_prob$prob),as.character(relevantProbesets))
我顺便检查了一下上面那个复杂的R函数跟我这次精简版的结果，发现这次才是对的，上面那个错了。
你们能发现上面那个为什么错了吗？

如果是在mysql里面它的表现形式如下；
mysql> select row_names,EGID,rowMeans from esetDataTable order by EGID limit 10;
+------------+-------+------------------+
| row_names | EGID | rowMeans |
+------------+-------+------------------+
| 38912_at | 10 | 85.8820246154773 |
| 41654_at | 100 | 301.720067595558 |
| 907_at | 100 | 273.100008206028 |
| 2053_at | 1000 | 354.207060199715 |
| 2054_g_at | 1000 | 33.8472900312781 |
| 40781_at | 10000 | 425.007640082848 |
| 35430_at | 10001 | 152.885791914329 |
| 35431_g_at | 10001 | 181.915087187117 |
| 35432_at | 10001 | 184.869017764782 |
| 31366_at | 10002 | 44.9716205901791 |
+------------+-------+------------------+
10 rows in set (0.05 sec)

如果要用SQL来做同样的事情需要下面这个语句，这个就非常简单啦！
select b.*
from test.esetDataTable b
inner join (
SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID
) as c on b.EGID=c.EGID and b.rowMeans=c.val
结果是：8640 rows in set (4.56 sec) 看起来mysql还没有R语言快速，但是这个mysql语句很明显也不是很高效，但是我的mysql水平有限。
稍微解释一下这个mysql语句，其中SELECT max(a.rowMeans) as val, a.EGID FROM test.esetDataTable a group by a.EGID类似于R里面的match_row=aggregate(dat,by=list(dat[,1]),max)，就是把每个entrez ID对应着的最大值提取出来成一个新的表
而inner join on b.EGID=c.EGID and b.rowMeans=c.val 就是我们的tmp_prob=merge(dat,match_row,by=c("EGID","rowMeans")) 根据两列来合并两个数据框

十 21

一个MIT的博士要离开学术圈，结果引发了上千人的热烈讨论（下）

Posted on by

Dr Petra EichelsdoerferFebruary 17, 2014 at 11:17 AM

My university lacked the funds to support me while applying for my first grant after my post-doc. So I was unable to re-write it after its first rejection. Although it was heartbreaking to walk away, I had a family to support. This was nearly four years ago. How many others have had to make the same decision I have? What will be the long-term consequences?

我们大学也是资金紧张，所以我直到博士后出站才拿到资金的第一笔科研基金。大约四年前把，当我第一次被基金会拒绝后，我实在难以鼓起勇气再次申请基金。尽管我必须去申请，因为我还得养家糊口。会有多少人曾与我经历过同样的处境呢，又有多少人与我一样做了同样的选择呢？ Continue reading →

十 20

一个MIT的博士要离开学术圈，结果引发了上千人的热烈讨论（中）

Posted on by

的确是上千人进行讨论，我这里分两次翻译网友们的讨论结果：

Lenny TeytelmanFebruary 15, 2014 at 7:53 AM

The reason I did not want to publish this - a single voice is invariably dismissed. So, I want to assemble in a central place as many essays like this from students, postdocs, and professors. The funding crisis will not be addressed until the severity of it is acknowledged and NIH, politicians, and scientists are alarmed enough. Please e-mail me your stories to lenny at zappylab dot com (whether new or published elsewhere). I will put together a site aggregating all of them.

以前我不发表这个观点的原因是-势单力孤。所以我希望等到足够多的类似的观点由不同的学生，博士后，甚至教授提出来之后，把它们整合在一起。学术圈的资金危机很难解决，除非NIH本身意识到问题的严重性，而且最好是政客，科学家们也对此有着足够的警醒。所以请把你的经历发邮件给我，不管你以前是否也在其它什么场合抱怨过，我将对他们做一个汇总，统一发表出来。

BioluminessaFebruary 17, 2014 at 10:55 AM

Hi Lenny. Thank you for sharing your perspective! Here is another to add to your collection of essays. I wrote this piece the other day coming to similar conclusions. http://bioluminate.blogspot.com/2014/02/the-seven-stages-of-grief-for-academic.html

谢谢你分享你的观点，这里有个故事可能比较符合你的要求。

sheiselsewhereFebruary 25, 2014 at 11:54 AM

I often get told that I shouldn't be so negative and that things will get better. Unfortunately, I don't have the time to wait. Here is my contribution.

http://sheiselsewhere.mosdave.com/2014/02/16/singing-for-supper/

长久以来，人们就告诉我我不应该如此的悲观，一切都会好起来的。不幸的是，我已经没有足够的时间去等待事情好转了。下面是我的经历。

Dregev21February 28, 2014 at 2:13 PM

Wow, thank you for posting this! I have gone through a very similar situation and have also decided to quit pursuing this dream. I was a 4th year PhD student at the University of Florida (where I had already had to change labs since my first mentor moved to UAB) and my project was going nowhere fast. I also started seeing academia for what it has become; an industry of cheap labor and false hopes. But like you, I stayed in it for as long as I could because of my love for science, learning and teaching. I quit and got out with a MS degree this past November and I am very happy with my decision. I began working as a research coordinator at UF, making more money and like you felt liberated and free from the constant stress of graduate work and research. I believe most students come in to graduate scholl for the same reasons, but it has become so disheartening and scary, that it didn't seem worth it to me anymore. I think it is important for current students to know and understand that there are other things to do in life that are more fruitfull, less stressful and just as intelectually stimulating and rewarding. In any case, thank you for sharing!

哇，谢谢你的分享！我有着同样的经历，也准备放弃一直追寻的梦想了。我是佛罗里达大学的一名在读博士生，这是我博士生涯的第四年，而博士期间，我已经换过一次导师了，因为我的前导师去了UAB，所以我的博士课题也换了。就我的经历来看，我也认为学生圈现在充斥着廉价劳动力和虚假的期望。然而就像你一样，凭着对科学的一腔热血，我坚持了很长时间。直到去年十一月，我辍学了，就拿到了一个硕士学位，并且我非常开心我做了这样的决定。我的第一份工作是在UF做一个研究协调员，可以挣很多钱了，也终于从无止境的毕业研究课题压力中解放出来了，跟你一样的自由了。我觉得是时候让现在的学生知道在我们的生命中还有这非常多的更有意义，会收获更大，但却压力更小的事情可以做。

Travels with Moby March 1, 2014 at 12:58 PM

Good Luck to you. I made this choice for many similar reasons about 6 years ago. But I was a college professor at a small college. You have aptly described the scenario. Even without the added stress of the grant machine, the choices that we are forces to make that divorce ourselves from family, friends to pursue this academic dream are incredibly costly. What I did find after a year in industry, is that I was not alone, I met former academics in industry and elsewhere that have expressed the same concerns. I wish you the best, and you are not alone.

祝你好运！我也做了同样的决定，在六年前。但是我曾经是一个小学院的一名教授。你非常精准的描述了我们这样的教授的状态。即使没有基金方面的压力，我们这样的选择也使得我们无法兼顾家庭，朋友。追求学术的道理牺牲真的好大。当我离开学术圈，进入产业界的第一年，我终于不再觉得孤单了，我遇到以前学术圈的好友的时候也听到过他们有着类似的抱怨。不管怎么样，还是希望你能好好的，毕竟，我们都不孤单。

Kevin ZelnioMarch 3, 2014 at 8:51 AM

Good luck! Life is better outside academia lol. I left 2 phds (got my masters after first one) and a decade long research career with 12 and then a 5+ year science communication career, left the country and started a microbrewery in Sweden. My skills as a scientist have been instrumental in my new profession as a beer maker (serious lab and sanitation skills here!) and a business person (improved and more diverse funding sources! AKA investors and people who drink BEER - which is like everybody). I cried a lot, I won't lie. Almost wrecked my marriage and the stress turned me into a horrible father for a while. Its just not a sustainable career for some types of people. Which is a shame, because the career is selecting for the same type of people and missing out on a diversity of life styles which could most likely benefit the scientific community in a number of ways. Here was my story: http://deepseanews.com/2013/02/19294/

祝你好运！学术圈外面的世界更精彩，O(∩_∩)O哈哈哈~我有着长达12年的科学研究经历和多于5年的科学普及经历，终于离开了我的祖国，来到了瑞典开始酿酒。在我还是科学家的时候掌握的各种技能对我现在作为一个啤酒酿造师还是蛮有帮助的，尤其是严肃认真的实验习惯以及无菌操作技巧。那时候申请各种稀奇古怪的基金所锻炼出来的沟通技巧也特别适合自己现在转行做商人。我曾经为这样的转变哭泣过，我承认。这几乎摧毁了我的婚姻生活，也在一段时间让我变成了一个糟糕的父亲。对某些人来说，这样的职业不是一辈子的事情，而且是某种意义上的丢脸。以下是我的故事。

UnknownMarch 3, 2014 at 2:10 PM

I don't see why I should view your departure as
a bad sign for the life sciences. As an engineer,
we celebrate when our students graduate, go
start a company or join an existing one, and
create products that make the world a better
place. Or, go work at a national laboratory,
the FCC, a non-profit, or any of the other types of
jobs where engineers make a contribution.

我不认为，应该把你的离开看作是科学界的损失。作为工程师，我们欢迎毕业了的同学开公司或者加入已有的公司，去创造一些让人们生活更美好的产品。或者，干脆去国家实验室工作，或者FCC也行，一个非盈利组织，或者可以参加很多其它类似的，对工程师有需求的工作

Lenny TeytelmanMarch 3, 2014 at 2:26 PM

First of all, it is terrific that you are supportive of graduate students who go on to be productive outside of academia! Unfortunately, in life sciences, you often lose support of your mentor the second you say that you do not plan to be a tenure track professor.
Second, and most importantly - the reason our departures and anxieties are cause for concern - being a professor, in the current funding climate, requires a level of sacrifice for science that fewer and fewer of the most talented and brightest scientists will make. Our taxpayers spend an extraordinary amount funding research. If the best scientists leave academia, research will suffer. Training of the future scientists will suffer. Science, inside and outside of academia will suffer.

首先，像你这样的毕业生，走出学术圈之后还能创造如此多的产品，必然是非常了不起的。不幸的是，在生命科学领域，你经常就会失去你导师的资助，就像你博文里面提到的第二点那样，你也不打算去争取一个终身教职了。

其次，最重要的是我们的离开仅仅是因为焦虑，在如今的基金资助条件下，成为一个教授，需要一定程度的牺牲，只有极少数的非常聪明，非常有才智的科学家能做到了。我们的纳税人还是投入了大笔资金在科学研究的，如果大量的优秀科学家都离开了学术圈，那谈何学术成果呢。一个科学家的训练是需要耗费大量时间和金钱的，他们都将因此受损，同时，不管是学术圈里面还是外面都会蒙受巨大的损失。

gregMarch 3, 2014 at 3:07 PM

Your story collection is a great idea. I hope you'll keep the sources of the site open? I bet a lot of people would like to contribute to making that project stand out - I would certainly be helping out.

Thanks a ton for your blog post. Your last point about leaving your wife if she'd treated you as badly as science does is awesome. I'm just coming to the realization that you seem to already have: the notion that "If you can see yourself possibly loving any profession as much as you love science, you're not cut out for science" is unhealthy - it's a mark of the sort of brainwashing that academia does to you.

Best wishes on your future path.

你这个想把所有类似的学术圈失意的故事收集起来的想法很棒。我希望你能持续收集，并且保持开放。我也确信会有非常多的人参与进来贡献他们的故事，让这个计划传播开来。我也会尽我的努力去推动它。

能看到你这个博文真的是三生有幸。你最后那句话（如果你妻子对你像科学对你那样你肯定把她给甩了）简直是太精彩了。我也开始有着你曾经有过的想法了，如果你对追求科学还抱有疑问，是因为你的爱不够坚定，这样的想法是不对的，这只是学术圈对你的洗脑。

Jessica WilsonMarch 6, 2014 at 11:43 AM

This is fantastic writing, despite the sadness. I sympathize (finishing PhD in neuroscience, considering heading out).

I'd love to try and make a video with some of the stories you've accumulated. I'm already looking through that Google Doc you posted right now, and my heart is breaking.

文章写的真棒，尽管让人感到莫名的悲伤和失望。作为一个刚刚结束了神经科学学习课程的博士生，我深有同感。

我很乐意为你收集的这些故事拍一个video，我也看完了你在google doc里面所表达的观点，实在是太震撼了。

CBMarch 12, 2014 at 6:09 PM

Really great stuff. I have reread this post a dozen times over the past couple weeks, as I am a postdoc currently on the precipice of throwing in the towel on my academic career. I find the last sentence particularly meaningful. I can't shake the feeling that giving up on this career that I have been laser-focused on for ten years feels an awful lot like a traumatic breakup. But the simple truth is exactly as you described, academic science simply doesn't respect its professionals nearly enough for the best of us to stick around.

Ugh, breakups suck!

 

Michael RuddyMarch 14, 2014 at 1:55 PM

How appropriate for a Valentine post ... if you do not love everything about what you are doing – move on until you find it!

O(∩_∩)O哈哈~在情人节发表这样的观点实在是太适合不过了。如果你实在是不喜欢你正在做的事情，果断的放弃，持续寻找直到找到你所爱的。

Nick EffordFebruary 15, 2014 at 8:53 AM

I sympathise and wish you a successful and fulfilling future, wherever that takes you. The pressures in UK academia are much the same, as is the relatively low pay. We've seen our pay fall 13% taking into account inflation over the last 5 or 6 years, and universities refuse to offer a decent pay increase despite increasing their income from students and despite the fact that they are sitting on huge cash reserves. My own institution would rather spend £50 million on new buildings than reward its staff for their dedication.

Like you and countless others, I'm reluctant to leave a job that can be very exciting and stimulating. But the truth is that the stress levels make it increasingly unsustainable. There is constant pressure to write papers and secure research funding and simultaneous pressure to improve teaching quality, but there is a failure to recognise that time is a finite resource, so one activity must inevitably be traded off against the other.

I don't expect to receive the same remuneration as I would in industry, but I do need one of two things to happen: either working conditions need to improve or the pay needs to improve to reflect the real pressures of the job. I've sacrificed too many evenings and weekends over the years, and that has had a negative impact on personal physical and mental health as well as family relationships. If something doesn't give soon, I could well end up following you out of academia.

The trade union for academics in the UK is currently locked in a bitter pay dispute with the universities. You can find out more about it at http://fairpay.web.ucu.org.uk/

Good luck!

Nick

我也深有同感，也期望你有个成功而且精彩的未来，不管你走向哪个领域。英国的学术圈压力也很类似，因为相对来说基金资助量都很小。在过去的5到6年间，我们的资金总量，考虑到通货膨胀，反而缩减了13%，大学却拒绝对研究基金做一个比较像样的提高，尽管他们向学生收取的学费更多了，而且他们的现金流也非常健康。我所在的研究所情愿花五千万欧元区修建一个大楼，也不会去奖励那些辛辛苦苦奉献着的教职工。

像你以及无数的其他类似经历的人一样，我也不情愿离开这个即兴奋又刺激的研究工作。但是压力的与日俱增让我的研究工作越来越难以为继。我们既要保证发表足够的文章，还要争取研究基金的支持，同时还要提高我们的教学质量，但是我们不得不承认，一个人的时间是有限的，所以我们必然会在有些方面做得不够。

我不期望能得到与工业界相当的报酬，我仅仅是期望我们的工作条件能得到改善，并且我们的报酬水平应该能对得住我们实际上所承受的工作压力。长年以来，我已经付出了无数个夜晚和周末，而这对我个人的身心健康是一个很大的影响，同时也极大的影响了我的家庭关系。如果这一切不尽快改善的话，我想，我应该马上就会追寻你的脚步，离开这学术圈了。

十 20

一个MIT的博士要离开学术圈，结果引发了上千人的热烈讨论（上）

Posted on by

看到了一篇好文，所以就抽空翻译了一下！比较搞笑的是这个博客在国内居然被墙了，所以我就只贴出我的译文吧。

再见吧，我的学术生涯！

在过去的12年间，我一直从事着科学研究和教学工作，而且我也乐于其中。然而，就在几个星期以前，我辞掉了MIT的博士后职位，放弃了长久以来追寻的学术梦想。随后，我觉得非常的自由快乐，而这对于美国的生命科学研究来说，的确是一非常不好的消息。

Michael Eisen，我的联合导师之一，他毕业于伯克利大学，最近在博客中写道：这是一个做科学研究非常好的时代，但对科学家来说确实一个非常烂的时代。几个月前，我与我的另一个联合导师Jasper Rine讨论了NIH研究基金会的资助风波。Jasper说道：除非NIH立刻醒悟，基金的管理方式需要重要的改变，否则，我们这一代的科学家命运将非常坎坷。而我就是这些命运坎坷的科学家的一员，所以我出局了。

2001年 我大学毕业之后，我拒绝了一个高薪的程序员职位。相反，我选择了在著名的冷泉港实验室做生物信息学，尽管薪水要低很多。我个人是非常兴奋能有这个机会用各 种计算技术来做生物学研究的。两年后，我如愿以偿的进入了该领域的研究生队伍，方向是分子生物学，与此同时，我的薪水在接下来的六年间都只有刚开始的一半 了。而在MIT做博士后的期间，我的薪水也没能回到十年前我作为一个初级程序员的时候，尽管那时的我技能欠佳，也没有什么拿得出手的专业领域知识。在商业领域中，个人报酬也当与他掌握的技能以及知识水平相当，但是，在学术领域，它们二者之间的关联度大打折扣。

幸运的是，我的妻子一直都很支持我对科学研究的热爱，从2006年起，她就在做物理系的助教来平衡我们的财政支持。她工作的报酬还行，所以我们能负担的起在剑桥的房租和日常开销。而月内我的另一个孩子即将来到这个世上，我还有一个孩子在托儿所，所以我必须得改变自己，来改善我们的财务状况，做一个居家好男人总比在MIT做一个日夜忙碌的博士后要强。

早在我还是 一名研究生的时候，我就意识到了追求学术生涯的种种弊端。我接受了微博的薪水，忍受了换学校单位的不稳定性以及教授们那近乎疯狂的工作量。我也接受了变幻 莫测的天气，在熬过来了十多年的研究生和博士后的日子，我本可以顺理成章的拿到教职。我甚至也能接受，在追寻着光荣而神圣的目标的过程中，五年后有可能会 被拒绝，从而不得不再次搬迁去另一个学校再寻找教职。我也能想象到，即使拿到了教职，我也不得不投入我的全部身心来为我的课题组争取研究基金。我能看到所 有的一切我将为我追求所爱而付出的代价。

然而，本来应该是五年之后选择作为一个教授的烦恼却提取让我害怕了，因为拿到NIH的资助的不确定性太大了。没什么好说的了，一想到我所提交的所有研究基金申请书都会被打道回府，我就恐惧万分。我本应该靠着研究基金来提升自己，但是现在却看起来更糟糕了。就我所认识的十多名在近两三年拿到教授职位的科学家，他们没有一个人拿到了NIH的巨大的资金支持，他们得到的只有拒绝，再拒绝，进一步的拒绝。其中有一个非常年轻，且富有才华的教授，申请基金已经高达13次了，但是也是一直被拒绝，尽管她的申请书写的非常棒，而且立题也非常新颖。因此，她马上就要失去她的实验室了，因为她的启动资金即将耗尽，而这一切对她来说都犹如噩梦般，甚至需要安眠药帮助度过。这一切，又如何能令我不害怕呢？

在青少年时 期，我就一直认为工作应该像度过周末一样愉快，而我也一直痴迷于此。在过去的十二年中，我一直试图在学术领域寻找这一点，而且认为只有学术这一条路能达到 这一点。幸运的是，一年前，我与一个生命科学家合作开发一个开放的，实时更新的中央资料库，我非常享受开发过程的每一个步骤，而且我也很确定，在公司也能 得到在科学界能得到的那种全心全意投入的感觉。

我现在还不 确定科学家们是否会用到我们所创造的产品，也不知道我在公司能否就有充足资金来追寻我的梦想。一个星期前的辞职，我的确是冒了很大的风险。风险是很大，但 这并不是疯狂的行为。真正疯狂的是按部就班的执着于学术圈。我可以这样说，通过这一年的创造各种科学性的产品，我比以前更接近于教授了。

我也明白， 很多读者都会认为我是一种吃不到葡萄就说葡萄酸的心态，或者认为我对学术的渴望并不是想象中的那么强烈，我真的企图是变得富有。如果这也是你所认为的，那 么你其实是抱着要科学家什么都不去想，只安心的做一个简简单单的教授的想法。我的确是热爱过我所从事的研究教学工作，我也好想念那些美好的日子，尽管想起 来很受伤。但我也爱我的妻子，如果她像学术界对待科学家那样对待我，我早就离开她了。

原文地址：http://anothersb.blogspot.com/2014/02/goodbye-academia.html

我又翻墙把原文拷贝出来了，如下：

Goodbye Academia

I have enjoyed research and teaching for the last twelve years. Yet, I have resigned from my postdoctoral position at MIT a week ago, giving up on the dream of an academic position. I feel liberated and happy, and this is a very bad sign for the future of life sciences in the United States.

Michael Eisen, my co-advisor from graduate school at Berkeley recently wrote that it is a great time to do science but a terrible time to be a scientist. A few months ago I was discussing with my other co-adviser Jasper Rine the crisis in NIH research funding awards (better known as "lottery"). Jasper said that unless NIH wakes up and there is a major restructuring, we will lose an entire generation of scientists. I am a member of this generation, and I am out.

In 2001, about to graduate from college, I turned down a programming position at a hedge fund. Instead, I chose to do bioinformatics at Cold Spring Harbor Laboratory for a much lower salary. I was excited about the possibilities of doing biological research using computational tools. Two years later, I enthusiastically entered graduate school in molecular biology, with my salary dropping by half for the next six years. As a postdoctoral researcher at MIT, I am not even back to earning what I did ten years ago as a junior programmer with no skills or domain-specific knowledge. In a commercial setting, my compensation would have kept pace with my knowledge and skills, but in academia, there seems to be a complete decoupling of the two.

Luckily, my wife has always been supportive of my passion for science and balanced my foolhardiness with a practical job as a physician’s assistant since 2006. She is well compensated, allowing us to pay off our loans and afford the monthly expenses in Cambridge. With a daughter in daycare and another child due in a month, we would certainly be in a better financial shape with me as a stay-at-home dad than a postdoctoral scientist at MIT.

Science has also meant wrenching moves across the country. In 2003, we moved to California for me to begin my graduate studies. We both love New York, and my wife was devastated to leave her family and friends. In 2009, after many tearful discussions, she agreed to move to Boston from California for my postdoc. The next move for a professor position would surely require moving to yet another new place in the country.

As a graduate student, I was well aware of all of the negatives of an academic career. I accepted the miniscule pay, the inability to choose where to live, and the insane workloads of professors. I accepted the uncertainty of whether, after 10-12 years as a graduate student and postdoc, I would actually get a job as a professor. I accepted that even after attaining this lofty goal, five years later, I could be denied tenure and would have to move to another university or go into industry. I accepted that even with tenure, I would have to worry my entire life about securing research funding for the lab. I saw all of these as the price to pay for doing something that I love.

However, one aspect of being a professor has been terrifying me for over five years now – the uncertainty of getting funding from NIH. No let me rephrase that. What is terrifying is the near-certainty that any grant I submit would be rejected. I have been waiting for the funding situation to improve, but it seems to only be getting worse. I personally know about ten scientists who have become professors in the last 3-4 years. Not a single one of them has been able to get a grant proposal funded; just rejection, after rejection, after rejection. One of these is a brilliant young professor who has applied for grants thirteen times and has been rejected consistently, despite glowing reviews and high marks for innovation. She is on the brink of losing her lab as her startup funds are running out and the prospect of this has literally led to sleepless nights and the need for sleeping pills. How can this not terrify me?

I have been obsessed since my teens with the idea that work should be something one desires to come back to after a weekend. For the last twelve years, being an academic was the only path I saw toward this. Fortunately, a year ago, I co-founded a startup to create an open, up-to-date, central protocol repository for life scientists. I have enjoyed every step of getting ZappyLab going, and I am certain that the company will give me the feeling that I still get from science - wanting to go into work every day.

I don’t know yet if scientists will use what we are building. I don’t know if we will be able to raise the capital needed to build what I dream of building. By resigning from my postdoc a week ago, I have done something very risky. Risky, but not crazy. What seems crazy is aiming to stay in the academic track. I say this despite having had the most scientifically productive year of my life; I am closer to getting a professorship than ever before.

I realize that many will dismiss my story as a tale of sour grapes, or say that my desire is not strong enough or my primary motivation is to get rich. If that is your position, you are simply hoping that future scientists will be unable to love anything other than being a professor. I do love research and teaching with every fiber of my being. I will miss them and it will hurt. But I also love my wife, and if she had treated me the way academia treats its scientists, I would have left her long ago.

他们对生物信息学的讨论有点类似于我国的科学网，但是我国的科学网博主水平很有限
十 16

最全面的转录组研究软件收集

Posted on by
能看到这个网站真的是一个意外，现在看来，还是外国人比较认真呀， 这份软件清单，能看出作者的确是花了大力气的，满满的都是诚意。from: https://en.wiki2.org/wiki/List_of_RNA-Seq_bioinformatics_tools
https://en.wiki2.org/wiki/List_of_RNA-Seq_bioinformatics_tools软件主要涵盖了转录组分析的以下18个方向，看我我才明白自己的水平的确没到家，印象中的转录组分析也就是差异表达，然后注释以下，最多分析一下融合基因，要不然就看看那些miRNA，和lncRNA咯，没想到里面的学问也大着呢，怪不得生物是一个大坑，来再多的学者也不怕，咱有的是研究方向给你。

    1 Quality control and pre-processing data
        1.1 Quality control and filtering data
        1.2 Detection of chimeric reads
        1.3 Errors Correction
        1.4 Pre-processing data
    2 Alignment Tools
        2.1 Short (Unspliced) aligners
        2.2 Spliced aligners
            2.2.1 Aligners based on known splice junctions (annotation-guided aligners)
            2.2.2 De novo Splice Aligners
                2.2.2.1 De novo Splice Aligners that also use annotation optionally
                2.2.2.2 Other Spliced Aligners
    3 Normalization, Quantitative analysis and Differential Expression
        3.1 Multi-tool solutions
    4 Workbench (analysis pipeline / integrated solutions)
        4.1 Commercial Solutions
        4.2 Open (free) Source Solutions
    5 Alternative Splicing Analysis
        5.1 General Tools
        5.2 Intron Retention Analysis
    6 Bias Correction
    7 Fusion genes/chimeras/translocation finders/structural variations
    8 Copy Number Variation identification
    9 RNA-Seq simulators
    10 Transcriptome assemblers
        10.1 Genome-Guided assemblers
        10.2 Genome-Independent (de novo) assemblers
            10.2.1 Assembly evaluation tools
    11 Co-expression networks
    12 miRNA prediction
    13 Visualization tools
    14 Functional, Network & Pathway Analysis Tools
    15 Further annotation tools for RNA-Seq data
    16 RNA-Seq Databases
    17 Webinars and Presentations
    18 References

 

十 16

3000多份水稻全基因组测序数据共享-主要是突变数据

Posted on by

感觉最近接触的生物信息学知识越多，越对大数据时代的到来更有同感了。现在的研究者，其实很多都可以自己在家里做了，大量的数据基本都是公开的， 但是一个人闭门造车成就真的有限，与他人交流的思想碰撞还是蛮重要的。

这里面列出了3000多份水稻全基因组测序数据，都共享在亚马逊云上面，是全基因组的双端测序数据，共3,024个水稻数据，比对到了五种不同的水稻参考基因组上面，而且主要是用GATK来找差异基因的。
而且，数据收集者还给出了一个snp calling的标准流程
我以前也是用这样的流程
SNP Pipeline Commands

1. Index the reference genome using bwa index

   /software/bwa-0.7.10/bwa index /reference/japonica/reference.fa

2. Align the paired reads to reference genome using bwa mem. 
   Note: Specify the number of threads or processes to use using the -t parameter. The possible number of threads depends on the machine where the command will run.

   /software/bwa-0.7.10/bwa mem -M -t 8 /reference/japonica/reference.fa /reads/filename_1.fq.gz /reads/filename_2.fq.gz > /output/filename.sam

3. Sort SAM file and output as BAM file

   java -Xmx8g -jar /software/picard-tools-1.119/SortSam.jar INPUT=/output/filename.sam OUTPUT=/output/filename.sorted.bam VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

4. Fix mate information

   java -Xmx8g -jar /software/picard-tools-1.119/FixMateInformation.jar INPUT=/output/filename.sorted.bam OUTPUT=/output/filename.fxmt.bam SO=coordinate VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE

5. Mark duplicate reads

   java -Xmx8g -jar /software/picard-tools-1.119/MarkDuplicates.jar INPUT=/output/filename.fxmt.bam OUTPUT=/output/filename.mkdup.bam METRICS_FILE=/output/filename.metrics VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

6. Add or replace read groups

   java -Xmx8g -jar /software/picard-tools-1.119/AddOrReplaceReadGroups.jar INPUT=/output/filename.mkdup.bam OUTPUT=/output/filename.addrep.bam RGID=readname PL=Illumina SM=readname CN=BGI VALIDATION_STRINGENCY=LENIENT SO=coordinate CREATE_INDEX=TRUE

7. Create index and dictionary for reference genome

   /software/samtools-1.0/samtools faidx /reference/japonica/reference.fa
   
   java -Xmx8g -jar /software/picard-tools-1.119/CreateSequenceDictionary.jar REFERENCE=/reference/japonica/reference.fa OUTPUT=/reference/reference.dict

8. Realign Target 

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T RealignerTargetCreator -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -o /output/filename.intervals -fixMisencodedQuals -nt 8

9. Indel Realigner

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T IndelRealigner -fixMisencodedQuals -I /output/filename.addrep.bam -R /reference/japonica/reference.fa -targetIntervals /output/filename.intervals -o /output/filename.realn.bam 

10. Merge individual BAM files if there are multiple read pairs per sample

   /software/samtools-1.0/samtools merge /output/filename.merged.bam /output/*.realn.bam

11. Call SNPs using Unified Genotyper

   java -Xmx8g -jar /software/GenomeAnalysisTK-3.2-2/GenomeAnalysisTK.jar -T UnifiedGenotyper -R /reference/japonica/reference.fa -I /output/filename.merged.bam -o filename.merged.vcf -glm BOTH -mbq 20 --genotyping_mode DISCOVERY -out_mode EMIT_ALL_SITES
十 16

NGS数据比对工具持续收集

Posted on by

无意中看到了这个网站，比wiki的还有全面和专业。搜集了现有还算比较出名的比对软件，并且列出来了，还做了简单评价，里面对比对工具的收集，主要是基于2012年的一个综述《Tools for mapping high-throughput sequencing data》，相信应该是有不少人都看过这篇综述的，其实生物信息初学者应该自己去文献数据库找点感兴趣的关键词的综述多看看，广泛涉猎总没有坏处的。

<img src="http://www.ebi.ac.uk/~nf/hts_mappers/mappers_timeline.jpeg" alt="Mappers Timeline" width="800">

Features Comparison

The following Table enables a comparison of mappers based on different characteristics. The table can be sorted by column (just click on the column name). The data was collected from different sources and in some cases was provided by the developers. For execution times and memory requirements we refer to the above mentioned review (supplementary data is available here).

The Data column indicates if the mapper is specifically tailored for DNA, RNA, miRNA, or bisulfite sequences.The Seq.Plat. column indicates if the mapper supports natively reads from a specific sequencing platform or not (N). The version column indicates the version of the mapper considered. Read length limits are showed in two columns: minimum read length (Min. RL) and maximum read length (Max. RL.). Unless otherwise stated the unit is base pairs. The support for mismatches and short indels is also presented including, when possible, the maximum number of allowed mismatches and indels: by default the value is presented in bases; in some cases the value is presented as a percentage of the read size; or as score, meaning that mapper uses a score function. The alignments reported column indicate the alignments reported when a read maps to multiple locations. The alignment column indicates if the reads are aligned end-to-end (Globally) or not (Locally). The Parallel column indicates if the mapper can be run in parallel and, if yes, how: using a shared-memory (SM) or/and a distributed memory (DM) computer. The QA (quality awareness) column indicates if the mapper uses read quality information during the mapping. The support for paired-end/mate-pair reads is indicated in the PE column. The Splicing column indicates, for the RNA mappers, if the detection of splice junctions is made de novo or/and through user provided libraries (Lib). The Index column indicates if the reads or/and the reference are indexed. The number of citations was obtained from Google Scholar on 13 June 2015.
十 16

根据染色体起始终止点坐标来获取碱基序列

Posted on by
这次要介绍一个非常实用的工具，很多时候，我们有一个染色体编号已经染色体起始终止为止，我们想知道这段序列是什么样的碱基。当然我们一般用去UCSC的genome browser里面去查询，而且可以得到非常多的信息，多到正常人根本就无法完全理解。但是我如果仅仅是想要一段序列呢？
诚然，我们可以下载3G的那个hg19.fa文件，然后写一个脚本去拿到序列，但是毕竟太麻烦，而且一般这种需求都是临时性的需要，我们当然想要一个非常简便的方法咯。
我这里介绍一个非常简单的方法，是基于perl的cgi编程，当然，不需要你编程了。人家UCSC已经写好了程序，你只需要把网页地址构造好即可，比如chr17:7676091,7676196 ，那么我只需要构造下面一个网页地址
hg38可以更换成hg19，dna?segment= 后面可以按照标准格式更换，既可以返回我们想要的序列了。
网页会返回 一个xml格式的信息，解析一下即可。
This XML file does not appear to have any style information associated with it. The document tree is shown below.
<DASDNA>
<SEQUENCE id="chr17" start="7676091" stop="7676196" version="1.00">
<DNA length="106">
aggggccaggagggggctggtgcaggggccgccggtgtaggagctgctgg tgcaggggccacggggggagcagcctctggcattctgggagcttcatctg gacctg
</DNA>
</SEQUENCE>
</DASDNA>
很明显里面的aggggccaggagggggctggtgcaggggccgccggtgtaggagctgctgg tgcaggggccacggggggagcagcctctggcattctgggagcttcatctg gacctg 就是我们想要的序列啦。
赶快去试一试吧
当然你不仅可以搜索DNA，还可以搜索很多其它的，你也不只是可以搜索人类的
See http://www.biodas.org for more info on DAS.
Try http://genome.ucsc.edu/cgi-bin/das/dsn for a list of databases.
X-DAS-Version: DAS/0.95
X-DAS-Status: 200
Content-Type:text
Access-Control-Allow-Origin: *
Access-Control-Expose-Headers: X-DAS-Version X-DAS-Status X-DAS-Capabilities

UCSC DAS Server.
See http://www.biodas.org for more info on DAS.
Try http://genome.ucsc.edu/cgi-bin/das/dsn for a list of databases.
See our DAS FAQ (http://genome.ucsc.edu/FAQ/FAQdownloads#download23)
for more information.  Alternatively, we also provide query capability
through our MySQL server; please see our FAQ for details
(http://genome.ucsc.edu/FAQ/FAQdownloads#download29).

Note that DAS is an inefficient protocol which does not support
all types of annotation in our database.  We recommend you
access the UCSC database by downloading the tab-separated files in
the downloads section (http://hgdownload.cse.ucsc.edu/downloads.html)
or by using the Table Browser (http://genome.ucsc.edu/cgi-bin/hgTables)
instead of DAS in most circumstances.

 

十 16

nature发表的统计学专题Statistics in biology

Posted on by

生物学里面，唯一还算有点技术含量，和有点门槛，就是生物统计了，而这也是绝大部分研究者的痛点，有能力的，可以看看nature上面关于统计学的专题讨论，而且主要是应用于自然科学的统计学讨论。

里面有几句统计学名言警句：
Statistics does not tell us whether we are right. It tells us the chances of being wrong.
统计学并不会告诉我们是否正确，而只是说明我们错误的可能性是多少。
Quality is often more important than quantity.
数据的质量远比数量要重要的多
The meaning of error bars is often misinterpreted, as is the statistical significance of their overlap.
Good experimental designs mitigate experimental error and the impact of factors not under study.
文章列表：
Research methods: Know when your numbers are significant
Scientific method: Statistical errors
Weak statistical standards implicated in scientific irreproducibility
The fickle P value generates irreproducible results
Vital statistics
Experimental biology: Sometimes Bayesian statistics are better
A call for transparent reporting to optimize the predictive value of preclinical research
Power failure: why small sample size undermines the reliability of neuroscience
Basic statistical analysis in genetic case-control studies
Erroneous analyses of interactions in neuroscience: a problem of significance
Analyzing 'omics data using hierarchical models
Advantages and pitfalls in the application of mixed-model association methods
Quality control and conduct of genome-wide association meta-analyses
Circular analysis in systems neuroscience: the dangers of double dipping
A solution to dependency: using multilevel analysis to accommodate nested data
How does multiple testing correction work?
What is Bayesian statistics?
What is a hidden Markov model?
下面的这些文章，其实就是我们正常课本里面统计学的知识点，但是放在nature杂志发表，就顿时高大上了好多
Points of significance: Importance of being uncertain
Points of Significance: Error bars
Points of significance: Significance, P values and t-tests
Points of significance: Power and sample size
Points of Significance: Visualizing samples with box plots
Points of significance: Comparing samples €”part I
Points of significance: Comparing samples part II
Points of significance:  Nonparametric tests
Points of significance: Designing comparative experiments
Points of significance: Analysis of variance and blocking
Points of Significance:  Replication
Points of Significance:  Nested designs
Points of Significance: Two-factor designs
Points of significance: Sources of variation
Points of Significance: Split plot design
Points of Significance: Bayes' theorem
Points of significance: Bayesian statistics
Points of Significance: Sampling distributions and the bootstrap
Points of Significance: Bayesian networks
A study with low statistical power has a reduced chance of detecting a true effect, but it is less well appreciated that low power also reduces the likelihood that a statistically significant result reflects a true effect. Here, we show that the average statistical power of studies in the neurosciences is very low. The consequences of this include overestimates of effect size and low reproducibility of results. There are also ethical dimensions to this problem, as unreliable research is inefficient and wasteful. Improving reproducibility in neuroscience is a key priority and requires attention to well-established but often ignored methodological principles.
十 16

生物信息学学者学习mysql之路

Posted on by
我一直都知道mysql其实很有用的， 哪怕是在bioinformatics领域。也断断续续的看过不少mysql教程，只是苦于没有机会应用。毕竟应用才是最好的学习方法，正好这些天需要用了，我就又梳理了一遍作为一个生物信息学学者，该如何学习mysql数据库。
然后再搜搜一堆技巧
差不多就可以开始啦。
我们不拿数据库来做网页，所以需要的仅仅是查询公共数据库的数据，当然，一般人都会选择直接去网页可视化的查询，或者去ftp批量下载后自己写脚本来查询，我以前也是这样想的，所以感觉mysql没什么用，因为它能做的， 我写一个脚本都能做到。但是任何事物能发展到如此流行的程度毕竟还是有它的优点的。
而在我看来，mysql的优点就是，不需要存储大量的文件信息，随查随用，如果我们想把数据库备份到本地，就要建立一大堆的文件夹，存放各种refgene信息呀，entrez gene信息呀，转录本，外显子等等各个文件夹，每个文件夹下面还有一堆文件，而且还要分物种存储，总之就是很麻烦，但是在数据库就不一样啦。
比如我们可以连接UCSC的数据库 （前提是你的机器里面可以允许mysql这个命令，而且你可以联网）
mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A
就这么简单， 你就用mysql远程登录了UCSC的数据库，可以show databases;或者use database hg19 ; 等等
里面有两百多个数据库，主要是多物种多版本，然后如果我们看hg19这个数据库，里面还有一万多个数据表，包含着hg19的全面信息。
还有很多其它的公共数据库可以练习
来自于：https://www.biostars.org/p/474/#9095

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29
ENSEMBL http://uswest.ensembl.org/info/data/mysql.html
GO http://www.geneontology.org/GO.database.shtml#mirrors

1000 Genomes: since June 16, 2011: http://www.1000genomes.org/public-ensembl-mysql-instance

mysql -h mysql-db.1000genomes.org -u anonymous -P 4272

Flybase has direct access to its postgres chado database.
http://flybase.org/forums/viewtopic.php?f=14&t=114
hostname: flybase.org port: 5432 username: flybase password: no password database name: flybase
e.g. psql -h flybase.org -U flybase flybase

mysql -h database.nencki-genomics.org -u public
mysql -h useastdb.ensembl.org -u anonymous -P 5306

你都可以登录进去看看里面有什么，也可以练习练习mysql的语法，但是增删改查种的查是可以用的
然后我们可以用R或者perl或者Python来连接数据库，也是蛮好用的， 我现在比较倾向于R
所以我就简单看了一下这个包的说明书，然后成功连接了
#Connect to the MySQL server using the command:
#mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A
#The -A flag is optional but is recommended for speed
library(RMySQL)
my.port="";
my.user="genome";
my.password="";
my.db="hg19";
#there are 203 databases,such as hg18,hg38,mm9,mm10,ce10
con <- dbConnect(MySQL(), host=my.host, user=my.user,dbname=my.db)
dbListTables(con) # there are 11016 tables in this hg19 database;

是不是很简单呀，只有你认真的学习，其实这些应用的东西都还是蛮简单的。

下面这本书也比较好，就讲了R或者perl或者Python来连接数据库，很全面

当然，如果想看mysql在bioinformatics方面的应用，下面还有很多学习资料
进阶版还可以看看具体事例，GO数据库的设计：http://geneontology.org/page/lead-database-schema
从这个来看，python要比perl 好很多http://www.personal.psu.edu/iua1/courses/files/2010/week15.pdf
十 16

居然还可以出售TCGA的数据，只有你稍微进行分析一下即可

Posted on by
亮瞎了我的双眼，原来还可以这样挣钱。
这个数据库的作者在2011年发了一篇如何寻找融合基因的文章：*Edgren, Henrik, et al. "Identification of fusion genes in breast cancer by paired-end RNA-sequencing." Genome Biol 12.1 (2011): R6.

然后基于此，把TCGA计划里面的所有癌症样本数据都处理了，并且得到了融合基因数据集，然后就以此出售

价格高达一万欧元，折合人民币七万多，一本万利，而且人家TCGA计划的数据的公开而且免费的，他做了二次处理就可以拿来挣钱，让我感觉很不爽。
到目前为止他们处理了TCGA计划里面的7652个癌症样本的数据，建立了一个囊括28种癌症的融合基因数据集，并且打包成了一个叫做FusionSCOUT 的产品来出售。
价格如下：

Pricing of FusionSCOUT datasets:

  • Single gene in one cancer set                        490€    /  580$ per dataset
  • Single gene fusions across all cancers          4900€  /  5800$ dataset
  • Individual cancer set                                       990 €   /  1250 $ per dataset
  • Full TCGA dataset                                          9900€  /  12500$ per dataset
该网站是这样介绍他们的产品的，号称有3500个研究团体已经使用了他们的数据，但是我感觉纯粹是吹牛，毕竟他这篇文献也就一百多的引用量，再说3500次购买，就这一个产品就能让他成为亿万富翁了，想想都觉得可怕。而且这网站这么烂，中国访问速度是渣渣，也就是相当于失去了中国的所有土豪客户了，怎么可能还有3500的销量，搞笑！

One of the latest therapeutics angles in the fight against cancer is fusion genes and their regulation. To aid in fusion gene research and reveal the multitude of gene fusion event in cancer samples MediSapiens has developed a proprietary FusionSCOUT pipeline for identifying fusion genes from RNA sequencing datasets.

Currently we have analysed 7625 tumour samples from the TCGA project building a fusion gene dataset covering 28 different cancers within the TCGA project which can be accessed through our FusionSCOUT product.

Using this pipeline, we have discovered 3930 samples with gene fusions with 9667 different fusion genes. We´ve discovered numerous novel gene fusions as well as new cancer types in which previously known fusions appear.

You can now purchase these gene fusions datasets with few mouse clicks and get the worlds most comprehensive gene fusions from cancer sets within days

FusionSCOUT cancer Reports

With FusionSCOUT you can access the full listings of all fusion genes in specific cancer datasets. Find new leads for possible cause of the cancer, examine the pathways that are affected by different fusions, stratify patients by shared fusion genes or search for potential target for drugs and companion diagnostics.

Once you purchase a FusionSCOUT dataset we will send you a detailed report with information on the fused genes, sample ID from the TCGA dataset, fusion frequencies across the dataset as well as fusion mRNA sequences and lists of protein domains present in the fusion transcripts.

By ordering the MediSapiens FusionSCOUT dataset, you´ll get:

  • A list of all gene fusions that involve your gene of interest, across all TCGA cancer types
  • TCGA sample ID: s of the for the samples with fusions
  • Exact exon junctions for the fusions, including alternatively spliced variants and data on whether reading frame is retained
  • Detailed list of protein domains retained in the fusion genes
  • cDNA sequence for the fusion mRNAs

Contact us to access the most up-to-date and comprehensive datasets of fusion gene events in different cancers!contact@medisapiens.com

Check out also our Fusion Gene Detection pipeline service for your samples!

Dataset missing? Email us and well add your favorite dataset to FusionSCOUT!

FusionSCOUT Cancer sets, March 2015

Cancer type Number of samples Number of fusion genes
Acute Myeloid Leukemia, LAML 153 69
Adrenocortical carcinoma, ACC 79 115
Bladder Urothelial Carcinoma, BLCA 273 473
Brain Lower Grade Glioma, LGG 467 309
Breast Invasive Carcinoma, BRCA 1029 3267
Cervical Squamous Cell Carcinoma and Endocervical Adenocarcinoma, CESC 195 190
Colon Adenocarcinoma, COAD 287 212
Glioblastoma multiforme, GBM 170 379
Head and Neck Squamous Cell Carcinoma, HNSC 412 386
Kidney Chromophobe, KICH 66 19
Kidney Renal Clear Cell Carcinoma, KIRC 523 217
Kidney Renal Papillary Cell Carcinoma, KIRP 226 145
Liver Hepatocellular Carcinoma, LIHC 198 317
Lung Adenocarcinoma, LUAD 456 991
Lung Squamous Cell Carcinoma, LUSC 482 1374
Lymphoid Neoplasm Diffuse Large B-cell Lymphoma, DLBC 28 18
Mesothelioma, MESO 36 26
Ovarian Serous Cystadenocarcinoma, OV 420 1166
Pancreatic Adenocarcinoma, PAAD 84 46
Pheochromocytoma and Paraganglioma, PCPG 184 83
Prostate Adenocarcinoma, PRAD 336 859
Rectum Adenocarcinoma, READ 85 74
Sarcoma, SARC 161 799
Skin Cutaneous Melanoma, SKCM 355 620
Stomach Adenocarcinoma, STAD 190 311
Thyroid Carcinoma, THCA 506 195
Uterine Carcinosarcoma, UCS 57 229
Uterine Corpus Endometrial Carcinoma, UCEC 167 422
十 14

几个国外出名的跟生物信息学相关的会议

Posted on by
会议列表如下：
ASGH会议-Annual Meeting of the American Society of Human Genetics
AGBT会议-Advances in Genome Biology & Technology (AGBT)
ASM会议-annual meeting of the American Society for Microbiology
ASHI会议-The American Society for Histocompatibility and Immunogenetics 
BOSC-生物信息开放会议：Bioinformatics Open Source Conference
ISMB/ECCB会议
ACMG会议-The ACMG Annual Clinical Genetics Meeting
annual Biology of Genomes (BoG) meeting at Cold Spring Harbor
以上排名不分先后，
一年一度的美国人类遗传学协会（ASHG）年会是遗传学界的盛事，也是目前规模最大的人类遗传学会议。2015年的年会于10月6-10日在马里兰州的巴尔的摩举行，吸引了6500多名科学家参与。他们将在会议上介绍和讨论人类遗传学各个方面的最新进展。
会议官网是：http://www.ashg.org/
非常隆重，也受到业界追捧！
会议ppt均可下载，但是要翻墙

https://storify.com/andrewsu/ashg14-speaker-slides

基因组生物学技术进展大会（AGBT）

中文介绍：

The 23rd Annual International Conference on Intelligent Systems for Molecular Biology (ISMB 2015)
14th Annual European Conference on Computational Biology (ECCB 2015)
这也是一个老牌会议了，会议官网是：https://www.iscb.org/
2015年的会议资料可以直接下载了：
ASM会议就比较水一点：
会议官网是：http://www.asm.org/

WHO WE ARE

The American Society for Microbiology (ASM) is the oldest and largest single life science membership organization in the world. Membership has grown from 59 scientists in 1899 to more than 39,000 members today, with more than one third located outside the United States. The members represent all aspects of the microbial sciences including microbiology educators.

The mission of ASM is to promote and advance the microbial sciences.

ASM accomplishes this mission through a variety of products, services and activities.

  • We provide a platform for sharing the latest scientific discoveries through our books, journals, meetings and conferences.
  • We help strengthen sustainable health systems around the world though our laboratory capacity building and global engagement programs.
  • We advance careers through our professional development programs and certifications.
  • We train and inspire the next generation of scientists through our outreach and educational programs.

ASM members have a passion for the microbial sciences, a desire to connect with their colleagues and a drive to be involved with the profession. Whether it is publishing in an ASM Journal, attending an ASM meeting or volunteering on one of the Society's many boards and committees.

Big parts of our everyday lives, from energy production, waste recycling, new sources of food, new drug development and infectious diseases to environmental problems and industrial processes-are studied in the microbial sciences.

Microbiology boasts some of the most illustrious names in the history of science--Pasteur, Koch, Fleming, Leeuwenhoek, Lister, Jenner and Salk--and some of the greatest achievements for mankind. Within the 20th century, a third of all Nobel Prizes in Physiology or Medicine have been awarded to microbiologists.

ASHI主要是免疫学相关的：

官网是：http://www.ashi-hla.org/

The ASHI 41st Annual Meeting site is now live, for the latest updates visit 2015.ashi-hla.org.
About ASHI

The American Society for Histocompatibility and Immunogenetics (ASHI) is a not-for-profit association of clinical and research professionals including immunologists, geneticists, molecular biologists, transplant physicians and surgeons, pathologists and technologists. As a professional society involved in histocompatibility, immunogenetics and transplantation, ASHI is dedicated to advancing the science and application of histocompatibility and immunogenetics; providing a forum for the exchange of information; and advocating the highest standards of laboratory testing in the interest of optimal patient care.

BOSC-生物信息开放会议

在wiki里面有详细的介绍：

The Bioinformatics Open Source Conference (BOSC) is an academic conference on open source programming in bioinformatics organised by the Open Bioinformatics Foundation. The conference has been held annually since 2000 and is run as a two-day satellite meeting preceding the Intelligent Systems for Molecular Biology (ISMB) conference.

annual Biology of Genomes (BoG) meeting

ACMG会议是临床相关的，报道的比较少
会议官网是；http://www.acmgmeeting.net/

ABOUT

The ACMG Annual Clinical Genetics Meeting provides genetics professionals with the opportunity to learn how genetics and genomics are being integrated into medical or clinical practice. The ACMG Annual Meeting Program Committee has developed a high caliber scientific program that will present the latest developments and research in clinical genetics and genomics

 

十 14

用wget批量下载需要认证的网页或者ftp站点里面的pdf文档

Posted on by
这是一个终极命令，比写什么爬虫简单多了
wget -c -r -np -k -L -p  -A.pdf --http-user=CS374-2011 --http-passwd=AlgorithmsInBiology http://ai.stanford.edu/~serafim/CS374_2011/papers/
我这里的例子是斯坦福大学的生物信息学算法课程里面推荐阅读的的所有pdf格式的paper
课程的网址是：http://ai.stanford.edu/~serafim/CS374_2011/
可以看到，这个网站推荐的文献分成8大类，本身这个网站打开就需要登录用户名和密码
--http-user=CS374-2011 --http-passwd=AlgorithmsInBiology
每一篇文献的单独地址是http://ai.stanford.edu/~serafim/CS374_2011/papers/Miscellaneous_topics/Self-assembly_of_DNA/self_healing_and_proofreading.pdf 类似的格式。
我这里简单解释一下这些参数的意思：
-c -r -np -k -L -p  -A.pdf
-c 断点续传
-r 递归下载，下载指定网页某一目录下（包括子目录）的所有文件
-nd 递归下载时不创建一层一层的目录，把所有的文件下载到当前目录（特殊要求会选择这个参数）
-np 递归下载时不搜索上层目录，如wget -c -r www.xxx.org/pub/path/
没有加参数-np，就会同时下载path的上一级目录pub下的其它文件 （所以一定要加上这个参数，不然会下载太多东西的）
-k 将绝对链接转为相对链接，下载整个站点后脱机浏览网页，最好加上这个参数
-L 递归时不进入其它主机，如wget -c -r www.xxx.org/
-p 下载网页所需的所有文件，如图片等
-A 指定要下载的文件样式列表，多个样式用逗号分隔
至于最后的--http-user=CS374-2011 --http-passwd=AlgorithmsInBiology 就是登录该课程网站需要的用户名和密码
十 12

生物信息学工程师在美帝的工资水平

Posted on by

今天逛论坛的时候，我看了一个宾夕法尼亚大学的生物信息学招聘启事：https://psu.jobs/job/60050

很有趣的是，我看到了他们的工资层级，而他们要招聘的生物信息学工程师的待遇是K，L级别的，也就是最低也是5万美金的年薪，折合成人民币还是蛮可观的，虽然我不是很清楚这个待遇在美帝属于什么样的水平，当然跟美帝的程序员肯定是没得比的，但是比国内的大部分程序员都还有好了。
Salary Band Minimum Midpoint Maximum
A $16,104 $23,748 $31,392
B $17,712 $26,124 $34,524
C $19,152 $28,728 $38,304
D $21,072 $31,620 $42,156
E $23,604 $35,400 $47,196
F $26,436 $39,660 $52,872
G $29,136 $44,412 $59,712
H $33,192 $50,616 $68,040
I $37,848 $57,696 $77,580
J $42,444 $65,772 $89,136
K $49,236 $76,308 $103,392
L $57,120 $88,524 $119,928
M $66,240 $102,672 $139,116
N $78,168 $121,152 $164,148
O $90,768 $142,968 $195,168
P $107,124 $168,696 $230,280
Q $126,396 $199,056 $271,728
R $151,668 $238,872 $326,088
A Bioinformatics Analyst position is available within the Bioinformatics Consulting Center at The Pennsylvania State University.
 The position is supported by the Huck Institutes for the Life Sciences and requires the candidate to work with multiple project investigators to design and implement computational pipelines for data produced by high throughput sequencing instruments and others, with particular emphasis on metagenomics and microbiome analyses.
 Responsibilities include the following: developing and/or maintaining existing software pipelines for analyzing high throughput sequencing data; identifying, evaluating and documenting new methodologies to support ongoing research needs; writing code and developing solutions to computational biology problems, with particular emphasis on microbiome and related samples. The Bioinformatics Analyst will become part of an interdisciplinary team composed of other bioinformatics staff, students and researchers and is expected to interact with other life scientists at Penn State and our international partner institutions in Africa and Asia to assist them with identifying research goals, analytical support needs, while carrying out computational data analysis as needed. It is anticipated that approximately 50% of your effort will initially be dedicated to providing bioinformatics support and microbiome analysis pipeline development for high-profile collaborative infectious disease surveillance research and training projects in Tanzania as well as other countries in East Africa and South Asia and may involve a limited amount of international travel (once per year). This job will be filled as a level 3, or level 4, depending upon the successful candidate's competencies, education, and experience. Typically requires a Master's degree or higher in a field of study with focus on computational research methods or higher plus four years of related experience, or an equivalent combination of education and experience for a level 3. Additional experience and/or education and competencies are required for higher level jobs. In-depth understanding of the computational analysis required for processing data from genomic technologies and their applications: Microbiome, metagenomics, RNA-Seq, genome assembly, genomic data visualization, or others. Expertise in handling and processing data in common bioinformatics formats; knowledge of available bioinformatics tools and genomic data repositories; proven track record of delivering bioinformatics solutions; demonstrated programming skills in one or more programming languages: Python, Perl, Java, C and/or numerical platforms: R, MATLAB, Mathematica. Experience handling large data sets generated from sequencing instruments. Excellent communication skills. This is a fixed-term appointment funded for one year from date of hire with excellent possibility of re-funding.