Category Archives: 生信组学技术

组学(Omics)主要包括基因组学（Genomics），转录组学(transcriptomics），全外显子组学(Whole-exome sequencing)，免疫组学(immune repetoire)，蛋白组学（Proteinomics），代谢组学（Metabolomics），免疫组学（Immunomics），糖组学（glycomics ）等

八 09

RNAseq数据完整生物信息分析流程第一讲之文献数据下载

Posted on by

我这里拿的是bioconductor里面最常用的airway数据，因为差异表达分析在bioconductor里面是重点，它们这些包在介绍自己的算法以及做示范的时候都用的这个数据。可以在GEO数据库里面看到信息描述：http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778  可以看到是Illumina HiSeq 2000 (Homo sapiens) ，75bp paired-end 这个信息很重要，决定了下载sra数据之后如何解压以及如何比对。也可以看到作者把所有的测序原始数据都上传到了SRA中心：http://www.ncbi.nlm.nih.gov/sra?term=SRP033351  ，这里可以在linux服务器上面写一个简单的脚本批量下载所有的测序数据，然后根据GEO里面描述的metadata把原始数据改名。

for ((i=508;i<=523;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP033/SRP033351/SRR1039$i/SRR1039$i.sra;done
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 $id;done

需要自己看SRA里面的数据记录，上面的脚本不难写出，然后因为是Illumina的双端测序，所以我们用fastq-dump --split-3命令来把sra格式数据转换为fastq，但是因为这里有16个测序数据，所以最好是同步改名，我这里用脚本批量生成改名脚本如下：

为了节省空间，我用了--gzip压缩，该文件名，用-A参数。

nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_untreated SRR1039508.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Dex SRR1039509.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb SRR1039510.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N61311_Alb_Dex SRR1039511.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_untreated SRR1039512.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Dex SRR1039513.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb SRR1039514.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N052611_Alb_Dex SRR1039515.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_untreated SRR1039516.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Dex SRR1039517.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb SRR1039518.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N080611_Alb_Dex SRR1039519.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_untreated SRR1039520.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Dex SRR1039521.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb SRR1039522.sra &
nohup ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --split-3 --gzip -A N061011_Alb_Dex SRR1039523.sra &

可以看到这里的16个样本来源于同样的4个人，是HASM细胞系，处理详情如下：

测序基础：
HASM细胞系-human airway smooth muscle，
The Illumina TruSeq assay was used to prepare 75bp paired-end libraries for HASM cells from four white male donors under four treatment conditions:
1) no treatment;
2) treatment with a β2-agonist (i.e. Albuterol, 1μM for 18h);
3) treatment with a glucocorticosteroid (i.e. Dexamethasone (Dex), 1μM for 18h);
4) simultaneous treatment with a β2-agonist and glucocorticoid
and the libraries were sequenced with an Illumina Hi-Seq 2000 instrument.
我们这里只是先根据fastq数据比对到参考基因组，然后计算每个样本的表达量即可，后续的分组计算差异表达，就需要个性化了。

下载的sra大小如下：

-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 04:21 SRR1039508.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 05:20 SRR1039509.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 06:14 SRR1039510.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 07:05 SRR1039511.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 08:07 SRR1039512.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 09:17 SRR1039513.sra
-rw-rw-r-- 1 jmzeng jmzeng 3.1G Aug 9 10:56 SRR1039514.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 11:56 SRR1039515.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 13:02 SRR1039516.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.6G Aug 9 14:16 SRR1039517.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.3G Aug 9 15:17 SRR1039518.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 16:05 SRR1039519.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 16:56 SRR1039520.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.4G Aug 9 17:57 SRR1039521.sra
-rw-rw-r-- 1 jmzeng jmzeng 2.0G Aug 9 18:46 SRR1039522.sra
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 19:28 SRR1039523.sra

解压后成双端测序的fastq数据如下：

 -rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.5G Aug 9 20:12 N052611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 20:44 N052611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 289M Aug 9 20:44 N052611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 951M Aug 9 20:59 N052611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 954M Aug 9 20:59 N052611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.7G Aug 9 20:53 N052611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.5G Aug 9 20:45 N061011_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:59 N061011_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 16M Aug 9 20:45 N061011_Alb.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.4G Aug 9 20:48 N061011_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 20:00 N061011_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 759M Aug 9 20:00 N061011_untreated.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.9G Aug 9 20:03 N080611_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 19:59 N080611_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 535M Aug 9 19:59 N080611_Alb_Dex.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 2.1G Aug 9 20:06 N080611_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.6G Aug 9 20:01 N080611_untreated_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_Alb_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:08 N61311_Alb_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.2G Aug 9 08:07 N61311_Dex_2.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_1.fastq.gz
-rw-rw-r-- 1 jmzeng jmzeng 1.3G Aug 9 08:09 N61311_untreated_2.fastq.gz

接下来所有的分析就基于此数据啦

 

 

八 07

对CHIP-seq数据call peaks应该选取unique比对的reads吗？

Posted on by

对于CHIP-seq数据处理完全是自学的，所以有很多细节得慢慢学习回来，这次记录的就是当我们把测序仪的fastq数据比对到参考基因组之后，应该对比对的结果文件做什么样的处理，然后去给peaks caller软件拿来call peaks呢？我看过博客 提到只保留比对质量值大于30的，也看过博客提到只保留unique比对的reads，我这里拿一篇公共数据测试了一下它们的区别！数据描述如下： Continue reading →

八 02

根据比对的bam文件来对peaks区域可视化

Posted on by

之前分析了好几个公共项目，拿到的peaks都很诡异，搞得我一直怀疑是不是自己分析错了。终于，功夫不负有心人，我分析了一个数据，它的peaks非常完美！！！可以证明，我的分析流程以及peaks绘图代码并没有错！数据来自于http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE74311，是关于H3K27ac_ChIP-Seq_LOUCY，组蛋白修饰的CHIP-seq数据，很容易就下载了作者上传的测序数据，然后跑了我的流程！https://github.com/jmzeng1314/NGS-pipeline/tree/master/CHIPseq Continue reading →

七 28

4种方式下载roadmap计划的所有数据

Posted on by
精选的129个细胞系，细胞系的介绍如下：http://www.broadinstitute.org/~anshul/projects/roadmap/metadata/EID_metadata.tab
对每个细胞系，都至少处理了5个核心组蛋白修饰数据，还有其它若干转录因子数据。
官网介绍的很详细，我就不翻译了：

Continue reading →

七 28

6种方式下载ENCODE计划的所有数据

Posted on by
DNA元件百科全书(Encyclopedia of DNA Elements, ENCODE)ENCODE计划的重要性我就不多说了，如果大家还不是很了解，可以直接跳到本文末尾去下载一下ENCODE教程，好好学习。该计划采用以下几种高通量测序技术来刻画了超过100种不同的细胞系或者组织内的全基因组范围内的基因调控元件信息。本来只是针对人类的，后来对mouse以及fly等模式生物也开始测这些数据并进行分析了， 叫做 modENCODE
  • chromatin structure (5C)
  • open chromatin (DNase-seq and FAIRE-seq)
  • histone modifications and DNA-binding of over 100 transcription factors (ChIP-seq)
  • RNA transcription (RNAseq and CAGE)

Continue reading →

七 26

用UCSC提供的Genome Browser工具来可视化customTrack

Posted on by
customTrack，我这里翻译为自定义的测序片段示踪文件，可以追踪我们的reads到底比对到了参加基因组的什么区域，或者追踪参考基因组的各个区域的覆盖度，测序深度！翻译自：http://genome.ucsc.edu/goldenPath/help/customTrack.html  这个非常有用！！！
UCSC提供的Genome Browser工具非常好用，可以很方便的浏览我们的测序数据在参考基因组的比对情况，由于定义好了一系列track的文件格式，用户可以非常方便的上传自己的track文件，但是如果用户超过48小时没有浏览自己的数据，UCSC会默认删除掉这些数据，除非用户已经保存在session里面。或者用户可以分享这些自定义的reads示踪文件customTrack。

Continue reading →

七 26

wig、bigWig和bedgraph文件详解

Posted on by
我们一般会熟悉sam/bam格式文件，就是把测序reads比对到参考基因组后的文件！bam或者bed格式的文件主要是为了追踪我们的reads到底比对到了参加基因组的什么区域，而UCSC规定的这几个文件格式(wig、bigWig和bedgraph)用处不一样，仅仅是为了追踪参考基因组的各个区域的覆盖度，测序深度！而且这些定义好的文件，可以无缝连接到UCSC的Genome Browser工具里面进行可视化！
这个网站提供了这几种数据格式的构造及转换脚本：http://barcwiki.wi.mit.edu/wiki/SOPs/coordinates
对SE数据，可以用macs2 pileup --extsize 200 -i $sample.bam -o $sample.bdg 把bam文件转换为bedgraph文件，不需要call peaks这一步骤。
而UCSC的ftp里面可以下载bedGraphToBigWig $sample.bdg ~/reference/genome/mm10/mm10.chrom.sizes $sample.bw 把bedgraph文件转换为bw文件，其余的转换工具都可以下载。

Continue reading →

七 14

ChIP-Seq文献数据重新分析解读第二例

Posted on by
这篇文章是朋友推荐的， 我觉得作为CHIP-seq学习材料再好不过了，所以推荐给大家。是全基因组范围的BRCA1和PALB2的转录共激活机制的探究。请务必先看我的CHIP-seq自学系列教程，跟着好好学习！数据如下：
GSM997540    BRCA1    SRR553473.sra    Read 18878514 spots
GSM997541    PALB2    SRR553474.sra    Read 17615498 spots
GSM997542    P_Ser2    SRR553475.sra    Read 35396009 spots

Continue reading →

七 13

ChIP-Seq文献数据重新分析解读第一例

Posted on by
文章是：Genome-wide maps of H3K4me2/3 in prostate cancer cell line LNCaP，数据在GEO可以下载。GSE20042，下面的所有分析，需要26G的空间。
作者想看看用 dihydrotestosterone （雄激素）处理了 cancer cell line LNCaP 这个细胞系之后，看看组蛋白甲基化修饰变化，主要是看H3K4me2和H3K4me3这两种组蛋白甲基化区别，分成三组，分别是处理前，处理后4H和16H，共有5个条件的数据，但是有7个fastq文件。
七 12

转录组 de novo流程–包括转录本完整注释

Posted on by

有网友咨询过对于没有参考基因组或者转录组的物种，如何做RNA-seq分析。我觉得这个问题太大了，而且我还真的对这个没有经验。但是我以前看到过一篇文献，里面提到过一个非常全面的转录组 de novo组装注释流程，所以我摘抄了文章里面的生物信息学处理部分，分享给大家： Continue reading →

三 15

融合基因检测软件-soapfusion

Posted on by

开发单位：华大，SOAP系列软件套装！

功能：检测合基因
优点：在现有的各种软件里面表现算是最好的
算法：是hash index，跟其它bwt算法不太一样
其它软件有： FusionSeq [21], deFuse [22], TopHat-Fusion [23], FusionHunter [24], SnowShoes-FTD [25], chimerascan [26] and FusionMap [27]
具体的算法我没看，因为只是有需求，正好有一些RNA-seq数据又想看看样本融合基因情况。所以就测试这个软件，通俗点说，融合基因原理其实很简单，如果有足够多的reads一部分比对到一个基因，另一部分比对到另一个基因，就可以说明它们两个基因发生了融合现象！如果是PE测序，那么更方便，左右两端reads比对情况也可以考虑。我就不多说废话了，直接上教程吧！
一，软件安装
下载压缩包，解压后即可使用！！！
推荐用最新版，然后看作者说明书的时候也要看清楚！
我反正好几次都搞糊涂了，最后联系了作者才搞明白，作者说他想更新到2.0版本，直接用HISAT的比对sam文件来做，但是还在筹备中，我觉得有点悬！
1
解压后是一堆perl程序，都在source目录下，source目录下面还有bin下面附带了几个第三方软件，包括bwa，blast和soap，最后都用得着！
有个很重要的问题，一定要软件自带的perl模块添加到perl的环境变量。不然那些perl程序运行会报错！
配置文件需要修改，就把几个目录放进去即可

二，输入数据准备

这里最重要的就是制作数据库！！！
作者给了非常详细的制作过程，我觉得还是不够清楚，所以再讲一遍！
首先下载5个文件：
6.5K Jun 15  2009 cytoBand.txt.gz
3.0G Oct 12  2012 hg19.fa
2.5M Mar 15 10:30 HGNC_Gene_Family_dataset
38M Feb  8  2014 Homo_sapiens.GRCh37.75.gtf.gz
202 Jan 19 16:07 HumanRef_refseg_symbols_relationship.list

文件下载地址，作者已经给出了！

我把这些文件都放在的当前文件夹下面的raw这个子文件夹，因为我要当前文件夹作为该软件的database文件夹！！！
然后运行命令！
我在SOAPfuse-v1.27文件下面运行：
perl ../SOAPfuse-v1.27/source/SOAPfuse-S00-Generate_SOAPfuse_database.pl  \
-wg raw/hg19.fa  -gtf raw/Homo_sapiens.GRCh37.75.gtf.gz  -cbd raw/cytoBand.txt.gz   -gf raw/HGNC_Gene_Family_dataset \
-rft raw/HumanRef_refseg_symbols_relationship.list \
 -sd ../SOAPfuse-v1.27 -dd ./

这一步耗时很长，4~6小时，创造了transcript.fa和gene.fa，然后还对他们建立bwa和soap的index，所以有点慢！

构建成功会有提示：
Congratulations!
You have constructed SOAPfuse database files successfully.
These database files are all stored in directory you supplied:
/home/jmzeng/biosoft/SOAPfuse/db_v1.27/
They are all generated based on public data files you supplied:
whole_genome_fasta_file:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/hg19.fa
gtf_annotation_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/Homo_sapiens.GRCh37.75.gtf.gz
Chr_Bandregion_file:       /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/cytoBand.txt.gz
HGNC_gene_family_file:     /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HGNC_Gene_Family_dataset
gtf_segname2refseg_list:   /home/jmzeng/biosoft/SOAPfuse/db_v1.27/raw/HumanRef_refseg_symbols_relationship.list
这些目录很重要，接下来制作配置文件会用得着！
To use these database files, just set the 'DB_db_dir' in config file as belowed:
DB_db_dir  =   /home/jmzeng/biosoft/SOAPfuse/db_v1.27
配置文件需要修改下面5个
DB_db_dir = /DATABASE_DIR/
PG_pg_dir = /TOOL_DIR/source/bin
PS_ps_dir = /TOOL_DIR/source
PD_all_out = /out_directory/
PA_all_fq_postfix = PostFix
其实你仔细阅读了说明书，你就知道该修改成什么样子了！
最后制作sample list文件
我这里只有一个sample,所以文件就一句话即可
test test test 100
所以我的有下面两个文件，都是为了顺应作者的需求我才搞了test/test/test这么无聊的东西！！！
/home/jmzeng/test_for_soapfuse/test/test/test_1.fq.gz
/home/jmzeng/test_for_soapfuse/test/test/test_2.fq.gz
如果你有多个sample需要一起运行，你就要仔细读作者的readme了，它把这个配置文件搞得特别复杂！！！

三，运行命令

如果文件都准备好了，运行命令非常简单！！
perl SOAPfuse-RUN.pl -c <config_file> -fd <WHOLE_SEQ-DATA_DIR> -l <sample_list> -o <out_directory> [Options]

运行的非常慢！！！

因为需要重新比对，知道

四，数据结果解读

结果，作者已经说的很清楚了，我就不多说了！
一 31

画基因的外显子覆盖度图

Posted on by
一般情况下，我们得到了测序reads在基因组的比对情况文件bam格式的，里面的信息非常多，如果我想特定的查看某个基因的情况，那么我们可以选择IGV等可视化工具，但它并不是万能的，因为即使是一个基因，它也会有多个转录本，多个外显子。
所以，我们可以画它的外显子覆盖图，如下：横坐标是外显子的长度，纵坐标是测序深度，每一个小图都是一个外显子
 DMD.NM_000109
根据这个图，我们就可以很明显的看出，DMD基因NM_000109转录本的1，10-17号外显子缺失，用IGV一个个的看这些外显子区域，是同样的结果！可能是芯片捕获不到，也可能是样本本身变异，造成的大片段缺失。但是这个图的信息就非常有用！
那么，我们该如何画这样的图呢？
首先，我们需要找到需要探究的基因的全部转录信息，及外显子信息！
在hg19_refGene.txt里面会有，在UCSC里面可以下载，新手可能会比较麻烦，实在不行你去annovar的目录也可以找到！
1
那么，我们根据这个信息，就可以判断该基因的起始终止位点啦
然后用samtools的depth命令去找这个基因的全部片段的测序深度信息
最后再格式化成下面的三列数据
第一列是该外显子的坐标，从1到该外显子的成都
第二列是该外显子在该坐标的测序深度，通过samtools的depth命令得到
最后一列是该外显子的标记，从exon:79一直倒推到exon:1，因为该基因在染色体的负链，所以外显子顺序是反着的！
1 84 exon:79
2 84 exon:79
3 84 exon:79
4 85 exon:79
5 85 exon:79
6 86 exon:79
7 85 exon:79
8 87 exon:79
9 89 exon:79
10 91 exon:79
11 92 exon:79
12 95 exon:79
13 96 exon:79
14 96 exon:79
15 99 exon:79
16 99 exon:79
17 97 exon:79
最后根据这个txt文档，用R语言，很容易就画出上面那样的图片了！
这里面的信息量还是蛮大的！
APOB.NM_000384

 

十一 01

WES（六）用annovar注释

Posted on by

使用annovar软件参考自：http://www.bio-info-trainee.com/?p=641

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample3.varscan.snp.vcf > Sample3.annovar

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample4.varscan.snp.vcf > Sample4.annovar

/home/jmzeng/bio-soft/annovar/convert2annovar.pl -format vcf4  Sample5.varscan.snp.vcf > Sample5.annovar

然后用下面这个脚本批量注释

image001

Reading gene annotation from /home/jmzeng/bio-soft/annovar/humandb/hg19_refGene.txt ... Done with 50914 transcripts (including 11516 without coding sequence annotation) for 26271 unique genes

最后查看结果可知，真正在外显子上面的突变并不多

23515 Sample3.anno.exonic_variant_function

23913 Sample4.anno.exonic_variant_function

24009 Sample5.anno.exonic_variant_function

annovar软件就是把我们得到的十万多个snp分类了，看看这些snp分别是基因的哪些位置，是否引起蛋白突变

downstream

exonic

exonic;splicing

intergenic

intronic

ncRNA_exonic

ncRNA_intronic

ncRNA_splicing

ncRNA_UTR3

ncRNA_UTR5

splicing

upstream

upstream;downstream

UTR3

UTR5

UTR5;UTR3

 

十一 01

WES（五）不同软件比较

Posted on by

主要是画韦恩图看看，参考：http://www.bio-info-trainee.com/?p=893

对合并而且过滤的高质量snp信息来看看四种不同的snp calling软件的差异

我们用R语言来画韦恩图

image001

可以看出不同软件的差异还是蛮大的，所以我只选四个软件的公共snp来进行分析

首先是sample3

image002

然后是sample4

image003

然后是sample5

image004

可以看出，不同的软件差异还是蛮大的，所以我重新比较了一下，这次只比较，它们不同的软件在exon位点上面的snp的差异，毕竟，我们这次是外显子测序，重点应该是外显子snp

image005

然后我们用同样的程序，画韦恩图，这次能明显看出来了，大部分的snp位点都至少有两到三个软件支持

所以，只有测序深度达到一定级别，用什么软件来做snp-calling其实影响并不大。

 

image006

 

十一 01

WES（四）不同个体的比较

Posted on by

3-4-5分别就是孩子、父亲、母亲

我对每个个体取他们的四种软件的公共snp来进行分析，并且只分析基因型，看看是否符合孟德尔遗传定律

image001

结果如下：

image002

粗略看起来好像很少不符合孟德尔遗传定律耶

然后我写了程序计算

image003

总共127138个可以计算的位点，共有18063个位点不符合孟德尔遗传定律，而且它们在染色体的分布情况如下

我检查了一下，不符合的原因，发现我把

chr1 100617887 C T:DP4=0,0,36,3 T:1/1:40 T:1/1:0,40:40 miss T:DP4=0,0,49,9 T:1/1:59 T:1/1:0,58:59 miss T:DP4=0,0,43,8 T:1/1:53 T:1/1:0,53:53 T:1/1:50

计算成了chr1 100617887 C 0/0 0/0 1/1 所以认为不符合，因为我认为只有四个软件都认为是snp的我才当作是snp的基因型，否则都是0/0

那么我就改写了程序，全部用gatk结果来计算。这次可以计算的snp有个176036，不符合的有20309，而且我看了不符合的snp的染色体分布，Y染色体有点异常

image004 image005

但是很失败，没什么发现！

 

 

十一 01

WES（三）snp-filter

Posted on by

其中freebayes,bcftools,gatk都是把所有的snp细节都call出来了，可以看到下面这些软件的结果有的高达一百多万个snp，而一般文献都说外显子组测序可鉴定约8万个变异！

image001

这样得到突变太多了，所以需要过滤。这里过滤的统一标准都是qual大于20，测序深度大于10。过滤之后的snp数量如下

image002

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample3.bcftools.vcf >Sample3.bcftools.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample4.bcftools.vcf >Sample4.bcftools.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if /INDEL/;/(DP4=.*?);/;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$1"}' Sample5.bcftools.vcf >Sample5.bcftools.vcf.filter

 

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample3.freebayes.vcf > Sample3.freebayes.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample4.freebayes.vcf > Sample4.freebayes.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next unless /TYPE=snp/;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]"}'  Sample5.freebayes.vcf > Sample5.freebayes.vcf.filter

 

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample3.gatk.UG.vcf  >Sample3.gatk.UG.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample4.gatk.UG.vcf  >Sample4.gatk.UG.vcf.filter

perl -alne '{next if $F[5]<20;/DP=(\d+)/;next if $1<10;next if length($F[3]) >1;next if length($F[4]) >1;@tmp=split/:/,$F[9];print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[1]:$tmp[2]"}'  Sample5.gatk.UG.vcf  >Sample5.gatk.UG.vcf.filter

 

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample3.varscan.snp.vcf >Sample3.varscan.snp.vcf.filter

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample4.varscan.snp.vcf >Sample4.varscan.snp.vcf.filter

perl -alne '{@tmp=split/:/,$F[9];next if $tmp[3]<10;print "$F[0]\t$F[1]\t$F[3]\t$F[4]:$tmp[0]:$tmp[3]"}' Sample5.varscan.snp.vcf >Sample5.varscan.snp.vcf.filter

这样不同工具产生的snp记录数就比较整齐了，我们先比较四种不同工具的call snp的情况，然后再比较三个人的区别。

然后写了一个程序把所有的snp合并起来比较

image003

得到了一个很有趣的表格，我放在excel里面看了看 ，主要是要看生物学意义，但是我的生物学知识好多都忘了，得重新学习了 image004

十一 01

WES（二）snp-calling

Posted on by

准备文件：下载必备的软件和参考基因组数据

1、软件

ps：还有samtools，freebayes和varscan软件，我以前下载过，这次就没有再弄了，但是下面会用到

2、参考基因组

3、参考 突变数据

第一步，下载数据

第二步，bwa比对

第三步，sam转为bam，并sort好

第四步，标记PCR重复，并去除

第五步，产生需要重排的坐标记录

第六步，根据重排记录文件把比对结果重新比对

第七步，把最终的bam文件转为mpileup文件

第八步，用bcftools 来call snp

第九步，用freebayes来call snp

第十步，用gatk     来call snp

第十一步，用varscan来call snp

下面的图片是按照顺序来的，我就不整理了

image001 image002 image003 image004 image005

image006 image007 image008 image009 image010 image011 image012 image013 image014 image015 image016 image017 image018 image019 image020 image021

 

九 25

RNA-seq流程需要进化啦！

Posted on by

Tophat 首次被发表已经是6年前

Cufflinks也是五年前的事情了

Star的比对速度是tophat的50倍，hisat更是star的1.2倍。

stringTie的组装速度是cufflinks的25倍，但是内存消耗却不到其一半。

Ballgown在差异分析方面比cuffdiff更高的特异性及准确性，且时间消耗不到cuffdiff的千分之一

Bowtie2+eXpress做质量控制优于tophat2+cufflinks和bowtie2+RSEM

Sailfish更是跳过了比对的步骤，直接进行kmer计数来做QC，特异性及准确性都还行，但是速度提高了25倍

kallisto同样不需要比对，速度比sailfish还要提高5倍！！！

参考：https://speakerdeck.com/stephenturner/rna-seq-qc-and-data-analysis-using-the-tuxedo-suite

九 24

用 GMAP/GSNAP软件进行RNA-seq的alignment

Posted on by

软件的解说ppt ：http://www.mi.fu-berlin.de/wiki/pub/ABI/CompMethodsWS11/MHuska_GSNAP.pdf

软件的下载地址： http://research-pub.gene.com/gmap/
有研究者认为这个软件的比对效果要比tophat要好，虽然现在已经多出来了非常多的RNA-seq的alignment软件，我还是简单看看这个软件吧，它本来是2005就出来的一个专门比对低通量的est序列，叫GMAP，后来进化成了GSNAP
step1：下载安装GMAP/GSNAP
是一个标准的linux源码程序，安装之前一定要看readme  ，http://research-pub.gene.com/gmap/src/README
解压进去，然后源码安装三部曲,首先 ./configu  然后make 最后make install
会默认安装在 /usr/local/bin 下面，这里需要修改，因为你可能没有 /usr/local/bin 权限,安装到自己的目录，然后把它添加到环境变量！
step2 ：准备数据
比对一般都只需要两个数据，一是索引好的参考基因组，另一个是需要比对的测序数据。
但是这个GSNAP，还需要对应的GTF注释文件。
首先需要参考基因组：虽然软件本身提供了一个hg19的参考基因组，并且已经索引好了Human genome, version hg19 (5.5 GB)(http://research-pub.gene.com/gmap/genomes/hg19.tar.gz) ，但是下载很慢，而且不是对所有版本的GSNAP都适用。所以我这里对我自己的参考基因组进行索引。
gmap_build -D ./ -d  my_hg19.fa
然后取ensemble下载hg19的gtf文件。
然后还需要把自己下载的gtf文件也构建索引，需要两个步骤
cat my_hg19.gtf |  ~/software/gmap-2011-10-16/util/gtf_splicesites > my_hg19.splicesites
cat  my_hg19.splicesites  |   iit_store -o my_hg19.gtf.index
然后拷贝需要比对的RNA-seq测序文件
step3: 运行程序
就是一步比对而已
gsnap
-D /home/jschnable/gsnap_indexes/
-d arabidopsisv10
--nthreads=50
-B 5
-s  /home/jschnable/gsnap_indexes/arabidopsisv10.iit
-n 2
-Q
--nofails
--format=sam temp.fastq
> results.sam
参数有点多，自己看看说明书吧http://qteller.com/RNAseq-analysis-recipe.pdf 讲的非常详细。