十一 07

mysql的table居然有最大列限制

想着把TCGA的RPKM值矩阵表格写入到mysql,然后做一个查询网页给生物学家,我下载的是所有TCGA收集的mRNA表达数据集数据集-GSE62944 ,共9264个癌症样本,和741个正常组织的表达数据。当我想写入癌症表达矩阵的时候,报错了:
Error in .local(conn, statement, ...) :
could not run statement: Too many columns
简单搜索了一下,发现是mysql有最大列的限制,但是我不是很懂计算机,所以没太看明白该如何调整参数使得mysql列限制扩充:http://dev.mysql.com/doc/refman/5.7/en/column-count-limit.html 所以就把癌症表达矩阵根据癌症拆分了,癌种数量如下:

Continue reading

十一 07

【直播】我的基因组(五):测试数据及参考基因组的准备

我的全基因组数据还没拿到,而且还会推迟,简单说(tu)明(cao)一下原因(还好当初为了避免广告嫌疑一直没说是哪个公司负责测序,反正用的是illumina的hiseqX10这个测序啦,所以可以尽情的吐槽)。

心烦意乱的吐槽线 Continue reading

十一 03

热图最佳实践-pheatmap

用pheatmap来绘图首先要安装这个包,它就一个功能,画出热图即可,号称是pretty heatmap,的确比其它的好用很多。我以前写过《一步一步学习heatmap.2》的教程,很简单的那种,所以就没有公布在博客上面,结果发现很多其它博客居然能先我一步发出。其实包括本次的pheatmap指南,都没什么好发,的在R里面也是傻瓜式出图,无法就是自己熟练一下参数而已,又不是开发一个包,没什么技术含量。我这里单独提一下pheatmap是因为它的确非常好用,将会是我画热图的不二之选。比如下面这个,是我最喜欢的: Continue reading

31

用samtools idxstats来对de novo的转录组数据计算表达量

de novo的转录组数据,比对的时候一般用的是自己组装好的trinity.fasta序列(挑选最长蛋白的转录本序列)来做参考,用bowtie2等工具直接将原始序列比对即可。所以比对 sam/bam文件本身就包含了参考序列的每一条转录本序列ID,直接对 sam/bam文件进行counts就知道每一个基因的表达量啦!

本来我是准备自己写脚本对sam文件进行counts就好,但是发现了samtools自带这样的工具:http://www.htslib.org/doc/samtools.html 

如果是针对基因组序列,那么这个功能用处不大,但是针对转录本序列,统计出来的就是我们想要的转录本表达量。 Continue reading

26

【直播】我的基因组(四):计算资源的准备

大家久等了,Jimmy的测序数据还没有拿到手!但是,工欲善其事必先利其器!所以jimmy在等待自己基因组的这段时间里,准备好了自己计算资源!鉴于会有不少同志们会跟着直播来自己动手分析一个公共的全基因组重测序数据(参考数据下载地址详见直播一),我在这里和大家分享一下计算机资源及软件的基本要求。 Continue reading

26

【直播】我的基因组(二):科研和临床分析调研

全基因组重测序,最大的分析点就是在于找到跟参考基因组不一样的地方(科研分析流程),然后通过各种公共数据库来注释这些不一样的地方(snv,indel,cnv,sv)(临床分析流程)。然而对于这些不一样的地方,就需要严格结合质量值、测序深度、重复性等因素进行进一步的过滤及筛选,过滤掉假阳性,从而进行下一步的分析。

那么我们一起来看看科研及临床一般都是进行怎样的分析吧! Continue reading

26

【直播】我的基因组(一):直播的目的及意义

直播的目的及意义
近年来随着“名人效应"的带动,基因测序也逐渐进入大家的视野。基因检测是在分子水平上对人体遗传密码进行破译,通过单核苷酸多态性和GWAS的分析对人体患病风险进行预测,从而进行预防干预及个体化治疗。随着二代测序技术的涌现,使得个人基因组测序成本逐年下降(2016年个人全基因组纯测序,30X,不到一万)。在这样的前景中,有一部分人开始考虑用基因检测来代替普通体检,从而预测预防各项疾病。最近有许多消息冒出来说测序体检报告如同天书看不懂,那么如何看懂测序呢?那让我们一起来看Jimmy的基因组直播。作为一名生信人,能亲自分析自己的测序数据可是再好不过了。Jimmy早在13年就产生了这样的想法,然而条件一直不是很成熟。但是执着的jimmy一直没有放弃这个念头,终于等到了天时地利人和的现在,决定开始和大家一起分享自己基因组的测序分析及解读的全过程。
Continue reading
22

使用trimmomatic对illumina数据做质控-去接头还有去除低质量碱基

因为一直拿到的是公司给的特别好的数据,所以没太关注质控这个问题,最近拿到了raw data,才发现其实里面的门道挺多的。前面都是用cutadapt这个python软件来去除接头的,但是它有一个弊端,需要自己指定接头文件。正好朋友推荐了trimmomatic,是java软件,所以直接Google找到其官网,然后下载二进制版本解压即可使用!
反正对我的illumina测序数据来说,直接用它就可以把raw data 变成 clean data啦!
1

Continue reading

15

R一大利器之对象的操作函数查询

对于生物出身的部分生物信息学工程师来说,很多计算机概念让人很头疼,尤其是计算机语言里面的高级对象。我以前学编程的时候,给我一个变量,一个数据,一个hash,我就心满意足了,可以解决大部分我数据处理问题,可事情远比想象之中复杂。因为很多高手喜欢用封装,代码复用,喜欢用高级对象。在R的bioconductor里面尤其是如此,经常会遇到各种包装好的S3,S4对象,看过说明书,倒是知道一些对象里面有什么,可以去如何处理那些对象,提取我们想要的信息,比如我就写过一系列的帖子:

Continue reading

15

用lumi包来处理illumina的bead系列表达芯片

表达芯片大家最熟悉的当然是affymetrix系列芯片啦,而且分析套路很简单,直接用R的affy包,就可以把cel文件经过RMA或者MAS5方法得到表达矩阵。illumina出厂的芯片略微有点不一样,它的原始数据有3个层级,一般拿到的是Processed data (示例), 当仍然需要一系列的统计学方法才能提取到表达矩阵。我比较喜欢用bioconductor,所以下面讲一讲如何用lumi包来处理这个芯片数据!

这个lumi包的使用代码和说明书都有,按部就班的学一遍就好了。
如果仅仅是分析数据,那么并不难,但是每个分析步骤后面都隐含着一系列的统计学方法,想彻底搞清楚他它们, 就很难了。

Continue reading

15

illumina的bead 系列表达芯片扫盲

表达芯片大家最熟悉的当然是affymetrix系列芯片啦,而且分析套路很简单,直接用R的affy包,就可以把cel文件经过RMA或者MAS5方法得到表达矩阵。illumina出厂的芯片略微有点不一样,它的原始数据有3个层级,一般拿到的是Processed data (示例), 当仍然需要一系列的统计学方法才能提取到表达矩阵。接下来我们首先讲一讲illumina的bead 系列表达芯片基础知识吧: Continue reading

14

转一个Python的安利文章咯!

来自于我们生信技能树论坛的超级版主bioinfo.dong的好文一篇,比较符合我博客的思想,就友情转发一下:

原文链接见:http://www.biotrainee.com/thread-379-1-1.html 

刚接触生信的同学大都有个困惑,知道生物信息可能需要编程,可是选择什么语言呢?有人会说perl啊,Python啊,R啊,java啊,等等等等。目的不一样,选择也不一样,你可以说语言都没有区别,达到目的就行,当然没问题。可是我们也要知道每种语言都有其独特优势,你可以用perl倒腾出矩阵运算,也可以画出想要的图,可是没有R专业;你也可以用R的正则表达式处理文本,可是perl或者Python做正则会更方便一些。这不是比较帖,只是从一个Python体验者的角度来说一下为什么选择Python。我目前的编程组合是Python+R+Shell Scripting。

这篇文章比较适合编程初学者,常年用perl的老司机们可以随便看一下,虽说perl和Python很像,有了一门的基础,学另一门就容易多了,可是真让一个用了几年perl的人彻底换Python还是比较困难的,主要还是习惯问题。最初做生信的人大都以perl作为常用脚本语言。我也是从perl开始的,当年为了申请出国读Bioinformatics,认真把小骆驼书看了一遍。来美国之后的第一个导师刚好是教perl的,我又跟着学了一次,看完导师推荐的《Unix and Perl to the Rescue》,算是巩固加第二次入门。之后一年基本都是用perl来处理数据。一个偶然的机会,同学说一起学学Python吧,听说很好用,于是就在网上找了个教程把题目刷了一遍。虽说入了门,可是每次项目赶时间的时候第一个想到的还是用perl来解决,所以入门很久也没啥长进,我亲爱的同学因为perl用的太好,虽然知道Python很好用,可始终没法狠心转过来,而我因为本身perl学得也只是半斤八两,纠结了一段时间也就彻底放弃perl了。

先说用了很长时间perl再用Python觉得不习惯的点。

(1)首先是动物园的书,《learning perl》真是入门的典范。再看《Learning Python》,几千页,那么厚,我到现在也没法认真看下去。
(2)另外perl语句比较简洁,几个符号就可以讲清楚的,Python可能需要几行,比如按行读取,perl只要while(<>)就可以,而最初学Python的时候,光这个问题就困扰了很久。再比如perl正则匹配的$1, Python是match.group(1)。perl的简洁伴随的缺点是可读性较差,自己的代码写完了都不想再看,更不要说别人写的。
(3)perl的正则表达式是真的非常厉害,我已经不记得是怎么厉害的了,就只记得Python的re module刚开始接触不太好用,不过现在已经感觉不出区别了。
(4)通常一个Python脚本需要很多modules,不熟悉之前会觉得很痛苦,perl就比较少用到,我总共也没用几次,一方面说明我的perl确实学得不好,另一方面可能也真是不太好用,看到就觉得麻烦。但Python的modules一旦熟悉了会大大提高工作效率。

重点说一下Python的优点。Python作为编程语言真正的优势比如面向对象编程(OOP),可移植/扩展/嵌入,强大的爬虫功能,APP开发,web开发等都不在讨论范围之内,只从最实用的角度做一下说明:
(1)简单,适合作为入门语言。很多时候觉得读Python的代码像是在读简单的英文,或者觉得pseudocode稍微一改就可以在Python里run了。Python还规范了很好的写作格式,该缩进的必须缩进,这样更增强了可读性。同时提高了代码重复利用的可能(很多时候perl代码写完就不想读了,三个月不用再回来已经看不太懂了,Python的就可以留着慢慢用。。。)
(2)Python社区活跃。有问题可以很容易搜索到解决方案。我perl的老师现在也转教Python了,问他为什么,他说perl的community不活跃,用Python是一种趋势
(3)作为开源语言,Python有很多非常好用的包,可以最大程度让我们避免把时间浪费在重复造轮子上。刚接触Python的时候我就觉得这简直是perl和R的整合,之前提过Python的scipy,numpy,pandas,matlibplot等等packages使其同样拥有了很强大的统计画图功能,我曾一度弃用R,用Python做所有的数据处理,数据分析和画图。不过现在又将这些工作交回了R,实验室本身是做统计的,用R显得入流一点:-)
(4)Python的jupyter notebook!!!这个是要强力推荐的!!!以前叫ipython notebook。用过R的都知道R Studio。jupyter notebook就是Python的Studio。以前写perl或者Python是不是这样的流程:写好了,存成.pl或.py格式,在shell里python xxx.py或者perl xxx.pl。运行完发现不好,有bug,打开文件找找bug在哪,再运行,还不行,唉,反反复复,好累。有了jupyter notebook你就可以边写边跑边改程序。有任何不确定的地方,都可以在notebook里直接测试,有任何bug都可以在notebook里直接改。简直方便到爆。现在用Anaconda安装jupyter还附赠很多包,方便又实惠。
(5)学好Python可以转行!!!跳出生物坑,奔向美好的互联网坑。前面提到的爬虫,APP开发,web编程都是很实用的技能。许多互联网公司也会专门招Python程序员,比如Google,比如Youtube,比如Dropbox。。。

我本专业是Bioinformatics,需要上一些计算机和统计的研究生课程,还记得算法课上老师第一节课就问,java和C++都会吧,如果不会的话Python总会吧,都不会的话这门课的作业写不了。就因为觉得自己还算会一点Python,把一次学习java的好机会浪费掉了

暂时就想到这么多。说的未必对。都是自己的体会吧。希望对初学者有用~

06

用cutadapt软件来对双端测序数据去除接头

一般来讲,我们对测序数据进行QC,就三个大的方向:Quality trimming, Adapter removal, Contaminant filtering,当我们是双端测序数据的时候,去除接头时,也会丢掉太短的reads,就容易导致左右两端测序文件reads数量不平衡,有一个比较好的软件能解决这个问题,我比较喜欢的是cutadapt软件的PE模式来去除接头!尤其是做基因组或者转录组de novo 组装的时候,尤其要去掉接头,去的干干净净!
cutadapt是经典的python软件,但是因为我的linux服务器有点问题 ,可能是root权限问题,没有用pip install cutadapt 安装成功,我懒得搞这些了,其实可以自己去下载cutadapt的源码,然后进入源码文件夹里面 python setup.py install --user 到自己的 ~/.local/bin下面。
所以我用conda安装了cutadapt软件,http://www.bio-info-trainee.com/1906.html 所以我需要 python ~/miniconda2/pkgs/cutadapt-1.10-py27_0/bin/cutadapt --help 才能调用这个软件,不过,问题不大,我也就是试用一下。 Continue reading

06

用sickle软件来对双端测序数据过滤低质量reads

一般来讲,我们对测序数据进行QC,就三个大的方向:Quality trimming, Adapter removal, Contaminant filtering,当我们是双端测序数据的时候,去除低质量的reads就容易导致左右两端测序文件不平衡,有一个比较好的软件能解决这个问题,而且软件使用非常简单! Continue reading

05

大型基因组拼装的乐高软件之MaSuRCA assembler使用指南

本文转载自 生信技能树 论坛特约作者MintMaSuRCA assembler 软件指导书,非常符合我博客的风格,也正式开启了我博客的转载之路。(前面的近400篇文章都是本人原创,手打,但是精力有限,以后文章更新频率会大大降低,但是会引入不少 技能树论坛特约作者的 好文!) Continue reading

05

用Miniconda,Bioconda来安装常见的生物信息学软件

这是生信技能树论坛的朋友推荐的,我试用了一下,的确非常方便~
生信最基础的基本功便是常用软件的安装和配置,但是不是所有软件都可以直接使用的[比如 annovar 等]。除了安装编译,有些软件所需环境的配置同样令人头疼。会不断报错提醒你那些东西没有安装。
    bioconda里面几乎涵盖了引用率较高的,好用的工具的打包资源,一键式安装,并且各自依赖的环境相互分隔,每次使用source activate env_name 来激活。使用source deactivate 来退出。具体软件列表见:https://anaconda.org/bioconda/repo 但是列表不支持搜索,可以去它的github里面去搜索  https://bioconda.github.io/
    而且不需要root权限也可以安装软件。

Continue reading

23

我挣大钱了?

最近跟一些志同道合的小伙伴们一起搭建了 生信技能树 的论坛,所以在社交上面加大了投入力度,认识了很多在生物信息学领域各个学习程度的同学,发现有些人问我的问题让哭笑不得,大意就是:我的生信菜鸟团博客有近四百篇文章,阅读10万+了,又是知乎上面的大V,现在又在建论坛,感觉生意很红火的样子,是不是挣了很多钱啊!

Continue reading

21

自学无参RNAseq数据分析第一讲之参考文献解读

这是我为新创办的 生信技能树 论坛写的帖子,也适合本博客,所以转载过来: http://www.biotrainee.com/thread-243-1-1.html 

以前做的都是有参转录组分析,只需要找到参考基因组和注释文件,然后走QC-->alignment-->counts->DEG-->annotation的流程图即可。
现在开始学习新的东西了,就是无参转录组分析,这里记录一下自己的学习笔记,首先还是资料收集,这次,我就针对性的看5个 全流程化的转录组 de novo 分析 文章,如下:
http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-554  2014年栀子花的花瓣衰老的标准de novo 转录组分析,数据如下:用Trinity做组装,用NCBI non-redundant (Nr) database库做注释,做了差异分析(栀子花花期分成4个阶段),GO/KEGG注释,然后做了RT-qPCR的实验验证。
多做了一个 Clusters of Orthologus Groups (COG)的数据库注释

Raw Reads
Clean Reads
Contigs
Unigenes
Annotated
Transcriptome
55,092,396
50,335,672
102,263
57,503
39,459

 

http://bmcgenomics.biomedcentral.com/articles/10.1186/1471-2164-15-236  2014 巴西橡胶树的研究,是一个综合多组织样本的RNA库,ployT建库,454测序,用的是est2Assembly 和gsassembler 软件做组装,用 NCBI RefSeq, Plant Protein Database 做注释,因为没有分组,所以不必做差异分析,只需要找SNV和SSR标记即可,最后也是做GO/KEGG注释

https://bmcgenomics.biomedcentral.com/articles/10.1186/s12864-016-2633-2 2015 萝卜,用illumina进行转录组测序,用Trinity组装,用RPKM值算unigene的表达量,也是用 BLASTx来对Trinity结果进行注释,注释到NR,NT,Swiss-Prot,GO,COG,kegg数据库,其中GO注释用的是Blast2GO,最后也做了RT-qPCR 实验验证,某些基因在leaf里面的表达量显著高于其它tissue,有原始数据:http://www.ncbi.nlm.nih.gov/sra/?term=SRX1671013
转录组分析结果结果:A total of 54.64 million clean reads and 111,167 contigs representing 53,642 unigenes were obtained from the radish leaf transcriptome.

http://www.nature.com/articles/srep08259 2015 芹菜 叶片发育中木质素的探究,测序的reads是A total of 32,477,416 quality reads were recorded for the leaves at Stage 1, 53,675,555 at Stage 2, and 27,158,566 at Stage 3, respectively.,也是用Trinity组装,kmer值设为25,组装结果:33,213 unigenes with an average length of 1,478 bp, a maximum length of 17,075 bp, and an N50 of 2,060 bp,然后用eggNOG/GO/KEGG数据库来注释。文章正文给了所用到的软件和数据库的详细链接
最后还用了 real-time PCR assays          来看 roots, stems, petioles, and leaf blade 这些组织的基因表达差异情况

http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0128659 对 三疣梭子蟹 的卵巢和睾丸的转录组研究,,也是标准的转录组de novo 分析流程,非常值得借鉴
NCBI有上传原始数据:SRR1920180  和SRR1920180  

总结好这5篇文献的数据分析流程,就差不多明白如何做无参的转录组de novo分析了