Daily Archives: 2015年3月23日

三 23

与hg19的突变相关的一些数据解释。

Posted on by

http://statgenpro.psychiatry.hku.hk/limx/kggseq/download/resources/

这个网站收集了大部分资料，我们就用它的，如果它倒闭了，大家再想办法去搜索吧。

与hg19的突变相关的一些数据解释210

其实这些文件都是基于NCBI以及UCSC和ensembl数据库的文件用一些脚本转换而来的，都是非常简单的脚本。

首先我们看看humandb/hg19_refGene.txt 这个文件，总共2.5万多个基因的共5万多个转录本。

与hg19的突变相关的一些数据解释325

19     可能是entrez ID，但是又不像。

NM_001291929    参考基因名

chr11   染色体

-

89057521

89223909

89059923

89223852

17      89057521,89069012,89070614,89073230,89075241,89088129,89106599,89133184,89133382,89135493,89155069,89165951,89173855,89177302,89182607,89184952,89223774,       89060044,89069113,89070683,89073339,89075361,89088211,89106660,89133247,89133547,89135710,89155150,89166024,89173883,89177400,89182692,89185063,89223909,

0

NOX4    基因的英文简称，通俗名

cmpl

cmpl

2,0,0,2,2,1,0,0,0,2,2,1,0,1,0,0,0,

然后我们看看hg19_snp141.txt这个文件

1       10229   A       -       .

1       10229   AACCCCTAACCCTAACCCTAAACCCTA     -       .

1       10231   C       A       .

1       10231   C       -       .

1       10234   C       T       .

1       10248   A       T       .

1       10250   A       C       .

1       10250   AC      -       .

1       10255   A       -       .

1       10257   A       C       .

1       10259   C       A       .

1       10291   C       T       .

1       10327   T       C       .

1       10329   ACCCCTAACCCTAACCCTAACCCT        -       .

1       10330   C       -       .

1       10390   C       -       .

1       10440   C       A       .

1       10440   C       -       .

1       10469   C       G       .

1       10492   C       T       .

1       10493   C       A       .

1       10519   G       C       .

1       10583   G       A       0.144169

1       10603   G       A       .

1       10611   C       G       0.0188246

1       10617   CGCCGTTGCAAAGGCGCGCCG   -

里面记录了以hg19为参考的所有的snp位点。

 

585

ENST00000518655 基因的ensembl ID号

chr1 + 11873 14409 14409 14409

4 基因有四个外显子

11873,12594,13402,13660, 12227,12721,13655,14409, 在基因的四个外显子的坐标

0

DDX11L1 基因的通俗英文名

none none -1,-1,-1,-1,

CTTGCCGTCAGCCTTTTCTTT·····gene的核苷酸序列

三 23

用annovar对snp进行注释

Posted on by

一、下载及安装软件

这个软件需要edu邮箱注册才能下载，可能是仅对科研高校开放吧。所以软件地址我就不列了。

用annovar对snp进行注释71

它其实是几个perl程序，比较重要的是这个人类的数据库，snp注释必须的。

用annovar对snp进行注释111

参考：http://annovar.readthedocs.org/en/latest/misc/accessory/

二，准备数据

既然是注释，那当然要有数据库啦！数据库倒是有下载地址

http://www.openbioinformatics.org/annovar/download/hg19_ALL.sites.2010_11.txt.gz

也可以用命令来下载

Perl ./annotate_variation.pl -downdb -buildver hg19 -webfrom annovar refGene humandb/

然后我们是对snp-calling流程跑出来的VCF文件进行注释，所以必须要有自己的VCF文件，VCF格式详解见本博客另一篇文章，或者搜索也行

http://vcftools.sourceforge.net/man_latest.html

三、运行的命令

首先把vcf格式文件，转换成空格分隔格式文件，自己写脚本也很好弄

perl convert2annovar.pl  -format vcf

/home/jmzeng/raw-reads/whole-exon/snp-calling/tmp1.vcf >annovar.input

用annovar对snp进行注释886

变成了空格分隔的文件

用annovar对snp进行注释900

然后把转换好的数据进行注释即可

./annotate_variation.pl -out ex1 -build hg19 example/ex1.avinput humandb/

 

四，输出文件解读

用annovar对snp进行注释1002

三 23

Snp-calling流程（BWA+SAMTOOLS+BCFTOOLS）

Posted on by

比对可以选择BWA或者bowtie，测序数据可以是单端也可以是双端，我这里简单讲一个，但是脚本都列出来了。而且我选择的是bowtie比对，然后单端数据。

首先进入hg19的目录，对它进行两个索引

samtools faidx hg19.fa

Bowtie2-build hg19.fa hg19

我这里随便从26G的测序数据里面选取了前1000行做了一个tmp.fa文件，进入tmp.fa这个文件的目录进行操作

Bowtie的使用方法详解见http://www.bio-info-trainee.com/?p=398

bowtie2 -x ../../../ref-database/hg19 -U  tmp1.fa -S tmp1.sam

samtools view -bS tmp1.sam > tmp1.bam

samtools sort tmp1.bam tmp1.sorted

samtools index tmp1.sorted.bam 

samtools mpileup -d 1000  -gSDf   ../../../ref-database/hg19.fa  tmp1.sorted.bam |bcftools view -cvNg -  >tmp1.vcf

然后就能看到我们产生的vcf变异格式文件啦！

 

当然，我们可能还需要对VCF文件进行再注释！

要看懂以上流程及命令，需要掌握BWA，bowtie，samtools，bcftools，

数据格式fasta，fastq，sam，vcf，pileup

 

如果是bwa把参考基因组索引化，然后aln得到后缀树，然后sampe对双端数据进行比对

首先bwa index 然后选择算法，进行索引。

然后aln脚本批量处理

==> bwa_aln.sh <==

while read id

do

echo $id

bwa aln hg19.fa $id >$id.sai

done <$1

然后sampe脚本批量处理

==> bwa_sampe.sh <==

while read id

do

echo $id

bwa sampe hg19.fa $id*sai $id*single >$id.sam

done <$1

然后是samtools的脚本

==> samtools.sh <==

while read id

do

echo $id

samtools view -bS $id.sam > $id.bam

samtools sort $id.bam $id.sorted

samtools index $id.sorted.bam

done <$1

然后是bcftools的脚本

==> bcftools.sh <==

while read id

do

echo $id

samtools mpileup -d 1000  -gSDf  ref.fa $id*sorted.bam |bcftools view -cvNg -  >$id.vcf

done <$1

 

 

==> mpileup.sh <==

while read id

do

echo $id

samtools mpileup -d 100000 -f hg19.fa $id*sorted.bam >$id.mpileup

done <$1