Category Archives: 未分类

四 02

使用bamdst完成公司外显子测序数据分析报告的重要环节

Posted on by

使用bamdst完成公司外显子测序数据分析报告的重要环节

目前主流的ngs科研服务，包括WES和RNA-seq价格都是透明的，反正建库四五百块钱，测序都是60元一个G左右，也就是说10G数据量的wes也就是一千块钱，同理RNA-seq也是如此，而且有趣的是标准的生物信息分析报告通常是附赠的，因为也的确是数据出来后必备的，只不过跟你的科研结论还差十万八千里。但是看懂这个报告又是必须滴，比如每个样本的下面这些指标 Continue reading →

四 02

坦白局

Posted on by

坦白局

春节旅途折腾，所以玩了个坦白局，很多人问的问题特别正能量，值得记录一下并且分享给所有的粉丝，因为是我语音回复，所以特别的口语化，而且我湖北口音比较重，首先感谢提问者帮忙把语音转换为文字版便于分享！ Continue reading →

四 02

StatQuest-2019

Posted on by

2018年有幸发现了一个非常棒的统计学学习视频，这个视频深入浅出的将复杂的统计学术语和理论解释的直白透彻，每一个统计学概念总是伴随着简单的不能再简单的图表和解释。

StatQuest是一个初学者极其友好的统计学教程

StatQuest作者总是不厌其烦的认真制作各个统计术语和理论的图表。 Continue reading →

四 02

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

Posted on by

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

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

四 02

有些时候警告会变成错误

Posted on by

有些时候警告会变成错误

经常就是打开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)
​
四 02

miRNA、LncRNA、CircRNA靠谱小结

Posted on by

miRNA、LncRNA、CircRNA靠谱小结

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

四 02

ANNOVAR 软件用法还可以更复杂

Posted on by

写在前面：

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

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

Continue reading →

四 02

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

Posted on by

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

首先需要理解：基因，转录本(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 →

四 02

cDNA 详解

Posted on by

cDNA 详解

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

四 02

MacOS的rjava加载失败问题

Posted on by

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 安装即可。

​