3大数据库超2万RNA-seq数据重新统一处理

各种大型计划产出的RNA-seq数据资源已经非常丰富了，但是大家都想把多个数据库联合起来分析，就不得不面对批次效应这个问题，所以UCSC团队就使用统一的流程把这些数据重新处理了，在亚马逊云上，一个样本花费1.3美元。
发表在：Nature Biotechnology publication: https://doi.org/10.1038/nbt.3772 Continue reading →

有些时候警告会变成错误

经常就是打开R或者终端写shell脚本，就会发现下面这样的警告：

During startup - Warning message:
Setting LC_CTYPE failed, using "C"
```<img class="wp-more-tag mce-wp-more" title="阅读更多…" src="data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7" alt="" data-wp-more="more" data-mce-resize="false" data-mce-placeholder="1" data-mce-src="data:image/gif;base64,R0lGODlhAQABAIAAAAAAAP///yH5BAEAAAAALAAAAAABAAEAAAIBRAA7">

通常我是选择无视它咯。

但也有时候，这个警告会成为致命错误。

比如安装一个GitHub来源的包：

```r
BiocManager::install(c('remotes','devtools'),ask = F,update = F) 
BiocManager::install('nik01010/dashboardthemes')

然后就报错了：

- preparing 'dashboardthemes':
v checking DESCRIPTION meta-information
- checking for LF line-endings in source and make files and shell scripts
- checking for empty or unneeded directories
- building 'dashboardthemes_1.0.2.tar.gz'

Error: (converted from warning) Setting LC_CTYPE failed, using "C"
Execution halted
Error in i.p(...) :
 (converted from warning) installation of package '/tmp/RtmpzlA16W/file6a8b5d45db7b/dashboardthemes_1.0.2.tar.gz' had non-zero exit status

这个时候大家会误以为是 had non-zero exit status 这个问题，所以大家求助就集中火力在这里，实际上是不对的。

很容易思考一下，就应该是去其它电脑运行同样的代码看看，

加上两句代码：

Sys.setenv(LANG="en_US.UTF-8")
Sys.setenv(LC_ALL="en_US.UTF-8")
BiocManager::install('ggpubr',ask = F,update = F)

就成功了：

- building dashboardthemes_1.0.2.tar.gz

* installing *source* package ‘dashboardthemes’ ...
** R
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded
* DONE (dashboardthemes)

miRNA、LncRNA、CircRNA靠谱小结

中心法则大家都不陌生，但是DNA <--> RNA —>蛋白质的传统的中心法则已经不能完整的概括基因遗传规律、解释所有的生命现象，随着非编码RNA的不断被发现 , 探测其功能已势在必行。相信在不久的将来，中心法则更完善的诠释生命的神奇。 Continue reading →

写在前面：

我多次在博客(生信菜鸟团)里面写过annovar软件的用法，而且在《直播我的基因组》里面也使用过它，https://mp.weixin.qq.com/s/GxqzKU7OPAMbjbZ5fPk7oQ 但那些都是最粗浅的使用而已，并没有深入了解ANNOVAR 软件的方方面面。

这次耗费15个小时系统性的回顾了该软件，希望可以做到教学上的最佳教程。

Continue reading →

Ubuntu16更新R的3.5版本

R到3.5因为引入了Bioconductor version: Release (3.8)，是一个破天荒地的改变，必须更新！

Ubuntu倒是很稳定，现在其实已经是Ubuntu18了。 Continue reading →

转录本融合位点上下游序列获取

首先需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列。

STAR-Fusion 结果解读

通常建议大家对RNA-seq数据使用 STAR-Fusion 来检测转录本融合现象，得到的结果如下：

文件 ‘star-fusion.fusion_predictions.abridged.tsv’, with the following format:

#FusionName JunctionReadCount SpanningFragCount SpliceType LeftGene LeftBreakpoint RightGene RightBreakpoint LargeAnchorSupport FFPM LeftBreakDinuc LeftBreakEntropy RightBreakDinuc RightBreakEntropy annots
THRA--AC090627.1 27 93 ONLY_REF_SPLICE THRA^ENSG00000126351.8 chr17:38243106:+ AC090627.1^ENSG00000235300.3 chr17:46371709:+ YES_LDAS 23875.8456 GT 1.8892 AG 1.9656 ["CCLE","FA_CancerSupp","INTRACHROMOSOMAL[chr17:8.12Mb]"]

可以看到列比较多，值得我们关心的是转录本融合的左右两个转录本的ID及其融合位点在基因组的坐标，如下：

THRA^ENSG00000126351.8 chr17:38243106:+ 
AC090627.1^ENSG00000235300.3 chr17:46371709:+

上面的坐标可以无限多，文件命名为 pos.txt 吧，后面写脚本需要用到，值得注意的是作者软件举例应该是hg19的基因组坐标，目前几乎都是hg38了，所以后面我脚本也是基于hg38的。

这个时候我们可能需要设计引物来对该融合转录本 进行验证，所以会需要这个融合点左右两个基因的指定转录本的cDNA序列。

我们很容易拿到各个转录本的基因组坐标，但是融合点的基因组坐标不能简单对应到转录本cDNA序列里面坐标。我们的突破点，就是找到融合点的基因组坐标到底对应到转录本cDNA序列的哪个位置。

物种的基因及其注释文件的下载

这里首选：https://asia.ensembl.org/info/data/ftp/index.html

wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/cds/Homo_sapiens.GRCh38.cds.all.fa.gz
wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.chr.gtf.gz

首先解析基因组注释文件，就是gtf的：

 1 havana exon 11869 12227 . + . gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; exon_number "1";

格式化代码如下：

zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -w havana |perl -F"\t" -alne '{next if $F[2] ne "exon";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > g2t2exon.txt
zcat Homo_sapiens.GRCh38.95.chr.gtf.gz |grep -w havana |perl -F"\t" -alne '{next if $F[2] ne "CDS";@a=split(/\"/,$F[8]); print join("\t",$a[1],$a[5],$a[9],$F[0],$F[3],$F[4],$F[4]-$F[3]+1) }' > g2t2cds.txt

结果如下：

ENSG00000223972 ENST00000456328 1 1 11869 12227 359
ENSG00000223972 ENST00000456328 2 1 12613 12721 109
ENSG00000223972 ENST00000456328 3 1 13221 14409 1189
ENSG00000223972 ENST00000450305 1 1 12010 12057 48
ENSG00000223972 ENST00000450305 2 1 12179 12227 49
ENSG00000223972 ENST00000450305 3 1 12613 12697 85
ENSG00000223972 ENST00000450305 4 1 12975 13052 78
ENSG00000223972 ENST00000450305 5 1 13221 13374 154
ENSG00000223972 ENST00000450305 6 1 13453 13670 218
ENSG00000227232 ENST00000488147 1 1 29534 29570 37

可以看到，一个基因会有多个转录本，然后每个转录本有多个外显子。不同的转录本的外显子有重叠。值得注意的是很多gtf里面记录的基因，在cDNA序列文件里面是不存在的，不过这个不影响我们的任务。

我上面两个代码分别得到的是外显子和CDS的坐标，后来发现CDS才是正确的， 因为并不是所有的外显子序列都会出现在cDNA序列里面。

首先基因组坐标转为转录本坐标

接下来需要写脚本把我们转录本融合位点那个基因组坐标，转为其转录本的相对坐标，这个时候普通的shell脚本已经无能为力，需要python或者perl这样的编程语言啦，就是把我们的 pos.txt 文件和自己制作的 g2t2exon.txt 关联起来，所以需要使用关联数组这个东西。

当然，也可以使用大家最擅长的R语言咯。

rm(list = ls())
options(stringsAsFactors = F)
a=read.table('input.txt')
head(a)
b=read.table('g2t2cds.txt')
colnames(b)=c('gene','transcript','exon','chr','start','end','exon_length')
head(b)
library(stringr)
re=NULL
tmp = apply(a,1,function(x){
 # x=a[1,]
 g=str_split(x[1],'[.]')[[1]][1]
 pos=as.numeric(str_split(x[2],':')[[1]][2])
 info=b[b[,1]==g,] ## one gene might has more than 1 transcripts.
 lapply(split(info,info[,2]), function(y){
 ## each transcript 
 apply(y,1,function(z){
 ## each exon. 
 if(z[5] <= pos & z[6] >= pos ){
 #print(c(z,pos))
 new = sum(as.numeric(y[y[,3] < as.numeric(z[3]),7])) + pos - as.numeric(z[5])+1 
 print(c(z,pos,new))
 re <<- rbind(re,c(z,pos,new))
 }
 }) ## each exon 
 }) ## each transcript 

}) ## each position 
colnames(re)=c('gene','transcript','exon','chr','start','end','exon_length','pos','new_pos')

代码非常复杂，需要一定编程水平才能理解。

1 ENSG00000121879 ENST00000643187 22 3 179234094 179235098 1005 179234680 3712
2 ENSG00000074800 ENST00000646370 13 1 8861007 8861429 423 8861330 1738
3 ENSG00000074800 ENST00000647408 12 1 8861002 8861429 428 8861330 1683
4 ENSG00000105976 ENST00000397752 1 7 116672392 116672577 186 116672464 73
5 ENSG00000105976 ENST00000456159 1 7 116672390 116672577 188 116672464 75
6 ENSG00000146648 ENST00000420316 7 7 55154011 55154152 142 55154029 1010
7 ENSG00000146648 ENST00000455089 6 7 55154011 55154152 142 55154029 888
8 ENSG00000077782 ENST00000526570 1 8 38427921 38430820 2900 38428753 833
9 ENSG00000037474 ENST00000502932 2 5 6603869 6604276 408 6604266 502

如上所述，融合位点的基因组坐标，成功转为了其对应的转录本序列坐标。

比如 ENST00000643187 这个转录本的坐标是

ENST00000643187.1 cds chromosome:GRCh38:3:179148574:179235098:1 gene:ENSG00000121879.4 gene_biotype:protein_coding

而融合位点在 179234680 ， 如果纯粹的使用它减去转录本起始坐标后是 86106 , 包含了大量的intron序列，所以需要找到其精准的外显子坐标。

比如这里的第22个外显子坐标是 3 179234094 179235098 ， 得到 586 的长度，再加上这个转录本前面的所有CDS的长度之和，最后是 3712 ， 就是该融合位点的转录本坐标啦。

然后根据转录本坐标及转录本序列获取

这个代码也不简单，需要读取我们下载好的cDNA序列文件：

rm(list = ls())
options(stringsAsFactors = F)
load(file = 're.Rdata')
fa=readLines('Homo_sapiens.GRCh38.cds.all.fa.gz')
fa=paste0( fa ,collapse = '')
fa=strsplit(fa,'>')[[1]]
fa=fa[-1]
tid=str_split(fa,'[.]',simplify = T)[,1]
seq=str_split(fa,']',simplify = T)[,2]
head(tid)
head(seq)
re=as.data.frame(re)
re=re[re$transcript %in% tid,]
re$tseq=seq[match(re$transcript,tid)]
head(re)
re$up=unlist(lapply(1:nrow(re),function(i){
 #i=1
 new_pos=as.numeric(re[i,9])
 l=nchar(re[i,10])
 if(new_pos>500){
 return(substring(re[i,10],new_pos-500,new_pos))
 }else{
 return(substring(re[i,10],0,new_pos))
 }

}))
re$down=unlist(lapply(1:nrow(re),function(i){
 #i=1
 new_pos=as.numeric(re[i,9])
 l=nchar(re[i,10])
 if((l-new_pos) >500){
 return(substring(re[i,10],new_pos,new_pos+500))
 }else{
 return(substring(re[i,10],new_pos,l))
 }

}))

判断前面转换好的转录本坐标是否扩充500后会越界，然后选取扩充500bp的序列即可。

值得注意的是转录本其实还有正负链信息是需要考虑的。

STAR-Fusion 提供工具组建融合后的转录本

其实并不需要自行写脚本进行探索，研究一下 STAR-Fusion ‘—denovo_reconstruct’ parameter. This requires that you include the ‘—FusionInspector inspect|validate’ setting. Continue reading →

cDNA 详解

我一直强调生物信息学工程师需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列。 Continue reading →

好课推荐：北京大学生科院《基因组学数据分析》

在李程老师的生物信息学平台社群看到了北京大学生科院《基因组学数据分析》，目录就勾起了我的学习欲望，独乐乐不如众乐乐，在得到李老师的同意后，现在在生信技能树平台分享这个课程。 Continue reading →

使用git删除文件

删除单个文件

如果只是删除本地的一个文件，通常是物理删除，然后git删除，再提交即可。 Continue reading →

MacOS的rjava加载失败问题

我的系统情况

RStudio Edition : Desktop
RStudio Version : v 1.1.383 (一直懒得更新)
OS Version : macOS 10.14 (18A391) (总是被迫更新)
R Version : 3.5.1 (2018-07-02)(现在可以更新到3.5.2了)

日常处理粉丝提问的时候，加载测试数据及代码

library(rJava)

很有趣的报错

> library(rJava)
arning: Error in : package or namespace load failed for ‘qdap’:
 .onLoad failed in loadNamespace() for 'rJava', details:
 call: dyn.load(file, DLLpath = DLLpath, ...)
 error: unable to load shared object '/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so':
 dlopen(/Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so, 6): Library not loaded: /Library/Java/JavaVirtualMachines/jdk-11.0.1.jdk/Contents/Home/lib/server/libjvm.dylib
 Referenced from: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rJava/libs/rJava.so
 Reason: image not found

解决bug思路

轻轻松松就谷歌搜索到前辈的经历：https://github.com/rstudio/rstudio/issues/2254

很简单，按照最高票答案，在我的mac终端输入：sudo R CMD javareconf 代码即可，但是二次报错

trying to compile and link a JNI program
detected JNI cpp flags : -I$(JAVA_HOME)/include -I$(JAVA_HOME)/include/darwin
detected JNI linker flags : -L/Users/jmzeng/Library/Java/Extensions -L/Library/Java/Extensions -L/Network/Library/Java/Extensions -L/System/Library/Java/Extensions -L/usr/lib/java -L. -ljvm
xcrun: error: invalid active developer path (/Library/Developer/CommandLineTools), missing xcrun at: /Library/Developer/CommandLineTools/usr/bin/xcrun
Unable to compile a JNI program

还是很明显，是xcrun问题，这个小毛病第十几次遇到了，每次更新macos系统都来一次，还是老规矩，使用代码；xcode-select --install 安装即可。

samtools和deeptools的numpy依赖冲突

本来我的conda环境里面安装的deeptools一直可以用，突然间就想安装samtools，发现之前的deeptools就开始报错 Continue reading →

基因长度之多少

无论RPKM或FPKM多么的遭人诟病，它的真实需求还是存在， 那么我们该如何合理的定义基因的长度呢？

这个时候需要理解：基因，转录本(transcripts,isoform，mRNA序列)、EXON区域，cDNA序列、UTR区域，ORF序列、CDS序列 这些概念，一个基因可以转录为多个转录本，真核生物里面每个转录本通常是由一个或者多个EXON组成，能翻译为蛋白的EXON区域是CDS区域，不能翻译的那些EXON的开头和结尾是UTR区域，翻译区域合起来是ORF序列，而转录本逆转录就是cDNA序列。 Continue reading →

经典的癌症外显子数据如何分析

癌症研究方法 :http://kidsgenomics.org/genomic-landscape-pediatric-cancers/ Continue reading →

感觉好像每隔几年学术界就掀起了一阵批判P值的讨论，正好我这里有个例子，分享给大家。
Continue reading →

因为TCGA计划跨时太长，这些年找somatic变异的软件也很多，所以TCGA团队下功夫在计划结束后（April 2018）完整的系统性的整理了最后的somatic突变数据。依托于文章：Scalable Open Science Approach for Mutation Calling of Tumor Exomes Using Multiple Genomic Pipelines March 201810.1016/j.cels.2018.03.002 Continue reading →

首先要了解什么是肿瘤异质性：http://cnvkit.readthedocs.io/en/latest/heterogeneity.html Continue reading →

你现在看到的是文献俱乐部2019年笔记
大名鼎鼎的TCGA计划回顾一下咯，关于AML研究发表在 N Engl J Med. 2013 May ， 算是很厉害了，附件一百多页的PDF详尽的描述了当时的数据分析方法，而且这个研究非常受重视，仅仅是关于 它的评论就有： Continue reading →

为了完成这个小探索，遇到了一个以前从来没有注意的问题，就是不同数据库对基因注释的记录差异问题。

• https://mp.weixin.qq.com/s/2QJvQxVECcxpJIsId1pHYA
• https://mp.weixin.qq.com/s/MDpX3tWQ7dojy84WjGuoZg Continue reading →

因为昨天看到了TxDb.Hsapiens.UCSC.hg38.knownGene 包来获取基因的坐标及长度跟其它主流数据库有差异，所以今天彻底比较一下TxDb.Hsapiens.UCSC.hg38.knownGene 包和CCDS数据库的差异。

详见：https://mp.weixin.qq.com/s/2QJvQxVECcxpJIsId1pHYA Continue reading →

本来以为做完了 http://www.bio-info-trainee.com/3727.html 准备工作，上课就问题不大，发现还是有一些很有意思的小问题，主要是计算机相关设置问题。 Continue reading →