十二 25

TCGA表达数据的多项应用之1–下载数据并且导入mysql

这个TCGA表达数据的多项应用系列帖子是应群里朋友的要求来写的,你们也可以继续提需求,我会接着写下去,其实从TCGA数据库里面下载到了数据之后,后面的所有分析都跟TCGA没有半毛钱关系了,大家要有这个想法,别三两句就问TCGA数据怎么分析,http://www.bio-info-trainee.com/?s=TCGA&submit=Search 本系列最后会形成一个shiny版本的交互式表达数据查询,处理,绘图,统计的网页APP。
我这里偷懒一下了,直接下载GEO里面的TCGA的表达数据,而不是去TCGA的官网里面下载:
它处理了目前(大概是2015年6月)TCGA收集的所有癌症样本的mRNA表达数据,并且统一处理成了count和RPKM两种表达量形式。 GEO地址:http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE62944

Continue reading

十二 19

【直播】我的基因组22:用IGV查看具体某个位点是否变异

下载IGV和导入文件的方法我就不多说了,可以直接在windows平台下使用,就跟你操作QQ一样,自己摸索就好了!

著名芬兰运动员Eero Mäntyranta,他拿过七枚奥运奖牌。他的血红细胞远超正常人水平,甚至一度被奥组委误以为服用了禁药。后来经过研究发现,他的EPOR基因上的一个位点rs121918116,发生了一个G>A的变异,使得他的血氧含量达到了普通人的150%,所以他耐力惊人。

在snpPedia里面可以查看这个位点的信息:http://snpedia.com/index.php/Rs121918116 Continue reading

十二 19

【直播】我的基因组(20):覆盖度详细探究

前面我们在第8讲提到了公司给我的一个报告的统计表格,有人反映不会做。

本来应该只需要给我6亿条reads的(PE150测序,人30X),但是足足给了我8.9亿条!(但事实上很多paper发表的基因组高于60X的也不少)

表格里面提到了好几个概念,比如duplicate的reads,一般说的PCR造成的duplicate,在找变异的时候需要去除掉。然后是那些比对到了不同染色体的reads pair,虽然只有2.29% ,也是需要重点分析的。(前面我也讲了如何提取以及分别分析它们!) Continue reading

十二 19

【直播】我的基因组(18):初步分析PCR duplication的情况

我博客里面有详细讲读原文查解去除PCR duplication的reads的原理和方法,还比较了samtools和picard这两个软件的区别,请点击阅看(仔细探究samtools的rmdup是如何行使去除PCR重复reads功能的),或者复制链接(http://www.bio-info-trainee.com/2003.html)到浏览器查看。 Continue reading

十二 19

【直播】我的基因组(17):初步分析一下multiple mapping 的情况

分析multiple mapping的情况,还是先用公司已经提供的bam文件来吧,后面我会用自己比对得到的bam文件再重新分析一次。

公司给我的bam数据是把有效测序数据通过 BWA(Li H et al.)比对到参考基因组 (b37),这是目前使用率很高的一个软件。在比对过程中,如果一个或一对 read(s) 在基因上可以有多个比对位置,BWA 的处理策略是从中选择一个最好的,如果有两个或以上最好的比对位置,则从中随机选择一个。这种多重比对(multiple hits)的处理对 SNP、INDEL 以及 CNV 等的检测有重要影响。 Continue reading

十二 19

【直播】我的基因组(16):提取左右端测序数据比对到不同染色体的PE reads

这类情况仅仅针对于双端测序数据,因为根据实验原理来看,对一个DNA片段,会把它的左右两端分别测序,但是测序仪器的测序长度有限,对本次实验来说,打断的DNA片段长度在350个碱基左右(这个长度只是一个分布,并不是真实值),理论来说测序是左右各150,加起来也就300,也就是说DNA片段中间还有50个碱基是测不到的(当然,实际上是有可能测通的)。而对这个配对的reads来说,来自于同一个DNA片段,所以理论上它们应该比对到同一条染色体的。也还是基于对sam格式的文件的理解,前面我们提到了sam文件的第3,7列指明了该reads比对到哪条染色体,以及该reads的配对reads比对到了哪条染色体(如果比对到同一条染色体,那么第7列是=符号)。所以我们只需要写脚本来提取即可! Continue reading

十二 19

【直播】我的基因组(15):提取未比对的测序数据

之前我们说了比对上的数据,那么会有人想到有没有没有比对上的数据呢?

既然是从我的血液里面提取到的DNA进行测序的,那么理论上测序仪出来的所有测序reads都应该是我的,也应该都可以比对到人类的参考基因组,但是实际过程中的确存在着未必对上的数据。

现有人类基因组毕竟只是个参考,也许我有某些独特的DNA序列呢?而且也不一定测序数据就都是人类的,也许我血液里面会有那么些微不纯粹呢?也有可能是某些片段变异的太多了,超过了常见的比对软件的承受能力,可以试用SHRIMP这个软件来提取一下。当然,这不是本讲的重点。 Continue reading

十二 19

【直播】我的基因组(14):bam文件给按照染色体给分割成小文件

昨天,我们了解了一下SAM格式的比对结果,不知道大家理解的怎么样。但是全基因组测序数据实在是太大了,即使比对后把sam文件压缩成二进制的bam文件也还有55G(如何压缩转换可查看直播十二),如果完整的导入IGV查看会略微考验计算机配置。

Continue reading

十二 19

【直播】我的基因组(13):了解sam格式比对结果

十一讲中将我们主要讲了如何将下机数据比对到参考基因组中。但是很多人对比对结果却是一头雾水。那我们现在来了解一下Sam格式的比对结果吧!

比对工具到现在已经多如牛毛了,见列表: https://en.wikipedia.org/wiki/List_of_sequence_alignment_software 。但是能被大多数人熟知的,就是bowtie和bwa(我们在十一讲中用的才是bwa),它们把测序数据比对到参考基因组之后,都会生成一个sam格式的文件。随后的大部分分析都是基于sam格式进行的分析,虽然Jimmy多次强调这些基础知识的重要性需要大家私下自学。但是由于这个sam文件实在是太重要了!!!所以,不得不亲自抽出一讲来说说它,后面也会基于此写十多篇文章: Continue reading

十二 15

制作自己的gene set文件给gsea软件

熟悉GSEA软件的都知道,它只需要GCT,CLS和GMT文件,其中GMT文件,GSEA的作者已经给出了一大堆!就是记录broad的Molecular Signatures Database (MSigDB) 已经收到了18026个geneset,但是我奇怪的是里面竟然没有包括cancer testis的gene set,MSigDB的确是多,但未必全,其实里面还有很多重复。而且有不少几乎没有意义的gene set。那我想做自己的gene set来用gsea软件做分析,就需要自己制造gmt格式的数据。因为即使下载了MSigDB的gene set,本质上就是gmt格式的数据而已:http://software.broadinstitute.org/cancer/software/gsea/wiki/index.php/Data_formats#GMT:_Gene_Matrix_Transposed_file_format_.28.2A.gmt.29 Continue reading

十二 15

CpG Islands记录文件下载的4种方式

这个也是读者来信最多的,关于基因组某些区域的起始终止坐标的下载问题,genomic feature的问题,一般是gtf文件或者bed文件,比如人类hg19上面的所有外显子的坐标记录文件,所有基因的坐标记录文件,所有lncRNA,rRNA等等,我这里拿CpG Islands记录文件下载的4种方式举例子给大家说明一下: Continue reading

十二 11

gene symbol 中的奇怪开头基因

这本是我为论坛的基础板块写的一个基础知识点,但是浏览量实在有限,不忍它蒙尘,特在博客重新发布一次!原帖见:http://www.biotrainee.com/thread-511-1-1.html

gene symbol 是非常官方的,由HUGO 组织负责维护,有专门的数据库HGNC database of human gene names | HUGO
以前分析数据的时候,有一些基因的symbol很奇怪,让我百思不得其解,比如
C orf 系列基因,
HS.系列基因,
KRTAP系列基因,
LOC系列基因,
MIR系列基因,
LINC系列基因
它们往往一个系列,就有好几百个基因;
C12orf44; Chromosome 12 Open Reading Frame 44;  这个是C orf系列基因的意思
MIR系列基因应该是 miRNA相关的基因
LINC系列基因应该就是long intergenic non-protein coding RNA
LOC系列基因,是非正式的,推定的,日后可能被更合适的名字替代
我这里做好了所有的基因对应关系,去生信菜鸟团QQ群里下载吧,共47938个基因的symbol和entrez gene id还有name,还有alias的对应!

1
还有一些RNA基因,根本就没有symbol,比如:CTA/B/C/D系列的
Aliases for ENSG00000271971 Gene
Quality Score for this RNA gene is 1
Aliases for ENSG00000271971 Gene
CTD-2006H14.2 5
External Ids for ENSG00000271971 Gene
Ensembl: ENSG00000271971
还有,如果你看到HS.开头的基因,它是unigene的ID了,已经不再是symbol啦。

十二 11

用R获取芯片探针与基因的对应关系三部曲-NCBI下载对应关系

这是系列文章,请先看:

用R获取芯片探针与基因的对应关系三部曲-bioconductor

ncbi现有的GPL已经过万了,但是bioconductor的芯片注释包不到一千,虽然bioconductor可以解决我们大部分的需要,比如affymetrix的95,133系列,深圳1.0st系列,HTA2.0系列,但是如果碰到比较生僻的芯片,bioconductor也不会刻意为之制作一个bioconductor的包,这时候就需要自行下载NCBI的GPL信息了,也可以通过R来解决:

##本质上是下载一个文件,读进R里面,然后解析行列式,得到芯片探针与基因的对应关系,看下面的代码,你就能理解了。 Continue reading

十二 10

解决阿里云博客的虚拟主机升级问题

首先感谢生信菜鸟团的各个小伙伴的大力支持,在阿里云的2年免费虚拟主机到期之后,我成功了续费了,但是坑爹的阿里云居然把我的IP地址和用户名都给替换了,导致了一些莫名其妙的bug。
虽然只是warning,不影响网站访问,但实在是影响界面美观,如下:
1

Continue reading

十二 09

【直播】我的基因组(十二):先粗略看看几个基因吧

昨天我们说到,测序得到的fastq文件map到基因组之后,我们通常会得到一个sam或者bam为扩展名的文件。SAM的全称是sequence alignment/map format,而BAM就是SAM的二进制文件。通常sam文件太大,我们会生成bam文件来节省空间。sam文件和bam文件的转换用samtools这个软件就可以完成。

samtools view -h abc.bam > abc.sam
samtools view -b -S abc.sam > abc.bam

我们已经拿到了bam文件,我这里就先用公司给我的bam文件吧,根据我的帖子:仅仅对感兴趣的基因call variation ,可以先了解几个比较有趣的基因的变异情况。我自己呢,对以下几个位点和基因比较感兴趣,就用他们来讲一下今天的内容吧!

1.STAT4上的rs7574865和HLA-DQ的 rs9275319是中国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因

2.V1aR基因是雄性标志性出轨基因。

3.GLI3和PAX1基因控制鼻孔的大小,而RUNX2基因控制鼻梁的宽度。DCHS2基因调控鼻子的突起程度,即决定鼻尖是否朝上和鼻尖的角度,或者说它决定了你的鼻子是否迷人挺拔。

4.肥胖有关的基因FTO(Fat Mass and Obesity Associated),最近发现了调控肥胖(主要是脂肪燃烧)的基因是IRX3 和IRX5。大约100个基因位点与BMI(身体质量指数)相关,600个基因位点与身高相关,160个基因位点与肥胖特征如腰臀比相关。6个新基因位点,这些位点位于LEMD2、CD47、GANAB、RPS6KA5/C14orf159、ANP32和ARL15基因内或周围。

那,我们就先关注这几个基因吧(不要问我为什么(-_-メ) )。

首先找到这些基因的坐标,看到如下:

其中V1aR基因这个雄性标志性出轨基因,在标准的基因命名系统里面其实是AVPR1A:http://www.genecards.org/cgi-bin/carddisp.pl?gene=AVPR1A ,这里面涉及到HUGO symbol的概念,这个genecard数据库也非常赞,基因相关信息都可以在这里面查找的。

有了这些坐标信息,我们就进入我们的基因组工作目录:

cd data/project/myGenome/

然后把坐标文件做好

因为公司给我的bam文件里面,用的参考基因组是GRCh37而不是hg19(两者区别在于chr是否标记),我们还是需要下载;

cd ~/reference

mkdir -p  genome/human_g1k_v37  && cd genome/human_g1k_v37

# http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/

nohup wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz  &

gunzip human_g1k_v37.fasta.gz

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.fai

wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/README.human_g1k_v37.fasta.txt

然后回到基因组工作目录,保证bam文件在上图中bamFiles那个目录,然后用下面这个脚本,批量提取我们感兴趣的基因的变异情况:

cat key_gene.list |while read id;

do

chr=$(echo $id |cut -d" " -f 1|sed 's/chr//' )

start=$(echo $id |cut -d" " -f 2 )

end=$(echo $id |cut -d" " -f 3 )

gene=$(echo $id |cut -d" " -f 4 )

echo $chr:$start-$end  $gene

samtools mpileup -r  $chr:$start-$end   -ugf ~/reference/genome/human_g1k_v37/human_g1k_v37.fasta bamFiles/P_jmzeng.final.bam  | \

bcftools call -vmO z -o $gene.vcf.gz

done

等三分钟就好了,结果如下:

前面我们说到有研究表明STAT4上的rs7574865和HLA-DQ的 rs9275319是国人群中乙型肝炎病毒(HBV)相关肝细胞癌(HCC)遗传易感基因,那么我们很容易去dbSNP数据库或者我最近强烈推荐 的snpedia数据库(吐血推荐snpedia数据库,非常丰富的snp信息记录)里面找到它的坐标。

6 32666295 :Rs9275319--HLA-DQ

2 191964633 :Rs7574865--STAT4

然后我检查了我刚才call到的variation文件,

zcat STAT4.vcf.gz |grep -w 191964633 显示为空。

zcat HLA-DQ* |grep 32666295  也是空。

哈哈,我完美的错过了这两个易感位点!!!!谢天谢地!!!

其余的我就不讲了,毕竟会涉及到隐私,我就讲这个方法吧!

文:Jimmy、吃瓜群众

图文编辑:吃瓜群众

十二 09

【直播】我的基因组(十):测序数据质量控制

质控之前我们在直播八的时候分析过,公司也给了我质控后的的数据,但是毕竟是别人做的,我们做为一个数据分析师,自己动手来验证一下公司给出的报告也是再好不过的了。大家可以跟着我先将下载数据进行一下质控。

因为此直播系列走得是半科普半技术路线,所以我这里show一个最常用也是最简单的测序质量控制软件,大名鼎鼎的fastqc软件,它是一个java软件,功能很单一,就是对你的测序数据生成一个网页版的可视化检测报告而已。这个软件的安装可以查看之前的直播贴(【直播】我的基因组(八):原始测序数据质量报告)。它在在linux或者windows平台都可以使用。直接下载这个压缩包: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.11.5.zip

我比较喜欢把这些软件放在biosoft文件夹下面(个人windows电脑),这个软件安装后会自带一些数据,大家感兴趣可以查看一下。

由于fastqc是免安装软件,直接解压后就可以直接使用。解压打开里面后缀是 .bat (相对于windows平台的批处理程序)的文件就打开fastqc啦,然后导入数据开始分析即可,静候一两个小时。

如果你用的是linux服务器,可以直接用unzip解压fastqc的zip压缩文件。里面有个fastqc的文件,就是fastqc的程序了。我们可以用fastqc  -o output dir [-(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN的命令让它进行质量控制。-o是用来指定输出文件的目录,注意是这里是不能自动新建目录的。输出的结果是.zip文件,默认自动解压缩,-noextract则不解压缩。-f用来强制指定输入文件格式,默认会自动检测。-c用来指定一个contaminant文件,fastqc会把overrepresented sequences往这个contaminant文件里搜索。后面加上你要质控的序列的文件名就可以了。

把所有的fastq.gz文件用fastqc软件处理得到的测序质量检测报告是一个html文件加上一个文件夹,如果没有解压缩需要用命令ls *zip|while read id;do unzip $id;done,把所有压缩包批量解压开。可以看到对每个测序数据它都进行了十几项统计结果和可视化的图片,对该款软件的结果感兴趣的可以下载(http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip) 文件,对原始数据处理前后的fastqc报告的区别显而易见。

然后批量抓取里面的%GC,Total sequences等信息,来跟之前公司给我的报告做比较,看看公司给我的结果有什么出入!

我以前写过帖子关于如何得到fastqc的统计表格:写脚本对fastqc的结果进行统计咯!


也就是说把多个qc的结果通过脚本整理在一起,方便查看。我们的统计结果如下:

当然一般不会有什么差别的,而且fastqc跑出来的结果都是合格的,公司对raw data得到clean的步骤仅仅是过滤掉不合格的reads,全部丢弃,而不是截断,豪气!!!

因为illumina的X10机器跑出来的数据一般都非常不错,我就没有在这里面下太多功夫,只是走个流程看一下测序质量,的确非常好,大家如果遇到质量比较差的数据,可以去我博客里面寻找各种解决方案。当然,质量控制不只是看序列的质量,还有很多小技巧,我会在后面的帖子里面专项讲解,比如我的数据是5条lane的数据合并起来的,那么lane的上样品是一定正确吗,那些没有比对上的reads是什么之类的相关问题。

请扫描以下二维码关注我们,获取直播系列的所有帖子!

1