搜索
查看: 3878|回复: 3

统计千人基因组计划里面的所有个体的纯杂合mutation比例

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-12-30 09:39:49 | 显示全部楼层 |阅读模式
在千人基因组计划的ftp服务器下载所有个体的突变信息文件如下:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502
  

它们这些文件本质上也是vcf格式的文件,只需要对vcf文件有足够的了解,就很容易从里面提取出纯合和杂合的信息,而且还区分出snp或indel,脚本略微有点复杂,我这里就不列出了。

重点就是要明白第9列后面的就是一个个样本,我们只需要关注 0|1     1|0     1|1 这3种情况就好了。
请自行阅读vcf说明书:VCF (Variant Call Format) version 4.1
GT : genotype, encoded as allele values separated by either of ”/” or “|”. The allele values are 0 for the reference allele (what is in the REF field), 1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on. For diploid calls examples could be 0/1, 1|0, or 1/2, etc.
但是值得注意的是”/” or “|”分割是不一样的,前者不需要care每个allele来自于哪条染色体,后者缺不一样。(应该是只有在大人群队列里面才能做到区分染色体来源的allele吧!)

脚本如下,我只给出一个例子:
[Perl] 纯文本查看 复制代码
 zcat ALL.chr10.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz |grep -v '^#' |head -10000 |perl -alne '{/VT=(.*?)\s/;$type=$1;if($type=~/SNP/){ $h{$_}{"SNP_$F[$_]"} +=1 foreach 9..$#F}else{ $h{$_}{"INDEL_$F[$_]"} +=1 foreach 9..$#F }}END{print $h{$_}{"SNP_1|1"}."\t".$h{$_}{"SNP_1|1"}."\t".$h{$_}{"SNP_1|0"}."\t".$h{$_}{"SNP_0|1"}."\t".$h{$_}{"SNP_0|0"}."\t".$h{$_}{"INDEL_1|1"}."\t".$h{$_}{"INDEL_1|0"}."\t".$h{$_}{"INDEL_0|1"}."\t".$h{$_}{"INDEL_0|0"}  foreach sort{$a <=> $b} keys %h}'

运算有点耗时,但是是准确无误的,唯一值得商榷的地方,就是两千多个个体在同一位点可能有不同的mutation情况,比如这个人是SNP,另一个人却是INDEL,所以这里的mutation type是多样性的, 我粗暴的用了匹配。
其实应该下载两千多人vcf文件分别来计算,而不是对这样的一个大文件统一处理,不过影响很小。

最后需要写成perl脚本对所有染色体的vcf文件批量统计,代码如下:
[Shell] 纯文本查看 复制代码
ls ALL*v5*gz  |while read id 
do 
echo $id
sample=`echo $id|cut -d"." -f 2` 
date
start=`date +%s`
zcat $id |grep -v '^#' |perl -alne '{/VT=(.*?)\s/;$type=$1;if($type=~/SNP/){ $h{$_}{"SNP_$F[$_]"} +=1 foreach 9..$#F}else{ $h{$_}{"INDEL_$F[$_]"} +=1 foreach 9..$#F }}END{print $h{$_}{"SNP_1|1"}."\t".$h{$_}{"SNP_1|1"}."\t".$h{$_}{"SNP_1|0"}."\t".$h{$_}{"SNP_0|1"}."\t".$h{$_}{"SNP_0|0"}."\t".$h{$_}{"INDEL_1|1"}."\t".$h{$_}{"INDEL_1|0"}."\t".$h{$_}{"INDEL_0|1"}."\t".$h{$_}{"INDEL_0|0"}  foreach sort{$a <=> $b} keys %h}'  > $sample.hom-het.mutation.stat
end=`date +%s`
date
runtime=$((end-start))
echo "Runtime for $sample was $runtime seconds"
done 


这个程序写的不好,耗时很夸张,希望有高手可以改善一下:ALL.chr10.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
Fri Dec 30 09:43:55 CST 2016
Fri Dec 30 11:27:42 CST 2016
Runtime for chr10 was 6227 seconds
ALL.chr11.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz
Fri Dec 30 11:27:42 CST 2016
Fri Dec 30 13:21:27 CST 2016
Runtime for chr11 was 6825 seconds
你看对不是那么大的10,11号染色体都耗费了两个多小时,运行完毕可是要两天呀!!!

结果如下:head chr10.hom-het.mutation.stat

72548        72548        55028        54497        3656187        11287        7919        7958        124688
73658        73658        54311        53699        3656602        11308        7735        7735        125136
70115        70115        54894        59717        3653514        10908        7797        8156        125029
73274        73274        50965        56620        3657457        11320        7671        8121        124746
73201        73201        53813        51206        3660108        11182        7759        7366        125607
72796        72796        53600        53510        3658420        11181        7723        7680        125352
70084        70084        52236        54699        3661313        10886        7390        7976        125666
71665        71665        54147        55254        3657247        11101        7619        8078        125058
74465        74465        54860        53772        3655198        11516        7898        7891        124597
73030        73030        53961        52981        3658328        11167        7876        7653        125167
可以看到对第10号染色体来说,第一个个体的SNP位点,1/1,1/0,0/1,0/0的个数分别是:72548        55028        54497        3656187(很明显我的代码有点疏忽了,把snp的1/1这种情况输出了两次)
其中0/0就不用考虑啊,就是没有snp位点在该个体。那么纯合和杂合的比例就很简单啦,72548/(55028+54497) 即可,每一行是一个个体,所有数据类似的处理,然后把所有的染色体综合起来考虑即可!
是不是很简单呀,这个数据就可以导入到R里面进行可视化。我把文件整理好了,在附件里面大家可以尝试着可视化一下:文件在附件里面,请下载。
然后R绘图代码是:
[AppleScript] 纯文本查看 复制代码
a <- read.table('merge.txt',stringsAsFactors = F)
a$snp_het_hom_ratio <- (a[,3]+a[,4])/a[,2]
a$nidel_het_hom_ratio <- (a[,7]+a[,8])/a[,6]
library(ggplot2)
dat <- a[,c(1,10,11)]
colnames(dat)=c('chr','snp','indel')
p=ggplot(dat,aes(x=chr,y=snp))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)

图片是:


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:ADDIS制作meta网状分析教程
下一篇:唯恩图 Venn或文氏图
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2017-1-7 15:13:25 | 显示全部楼层
本帖最后由 惠真-市二医院 于 2017-1-7 15:19 编辑

perl脚本!
我计算的百分率也包含了野生型个数:
[Perl] 纯文本查看 复制代码
  1 #!/home/huizhen/bin/perl
  2
  3 use strict;
  4 use warnings;
  5 use List::Util qw(sum);
  6
  7 die "perl $0 <vcf> <out>:" unless @ARGV == 2;
  8
  9 my $vcf = shift @ARGV;
 10 my $out = shift @ARGV;
 11
 12 open VCF,"<",$vcf or die "cant open:$!";
 13 open OUT,">>",$out or die "cant open:$!";
 14
 15 my $mut_num = 0;
 16 my $mut_other = 0;
 17 my $chr;
 18 my %hash;
 19 while(<VCF>){
 20     chmod;
 21     if(! /^#/){
 22         my @temp = split;
 23         $chr = $temp[0];
 24         if(/VT=SNP/){
 25             for(my$i=9;$i<=$#temp;$i++){
 26                 if($temp[$i] eq "0|0"){
 27                     $hash{"snp"}{"ref"} += 1;
 28                 }elsif($temp[$i] eq "0|1" or $temp[$i] eq "1|0"){
 29                     $hash{"snp"}{"het"} += 1;
 30                 }elsif($temp[$i] eq "1|1"){
 31                     $hash{"snp"}{"mut"} += 1;
 32                 }
 33             }
 34             $mut_num ++;
 35             }elsif(/VT=INDEL/){
 36                 for(my$i=9;$i<=$#temp;$i++){
 37                 if($temp[$i] eq "0|0"){
 38                     $hash{"indel"}{"ref"} += 1;
 39                 }elsif($temp[$i] eq "0|1" or $temp[$i] eq "1|0"){
 40                     $hash{"indel"}{"het"} += 1;
 41                 }elsif($temp[$i] eq "1|1"){
 42                     $hash{"indel"}{"mut"} += 1;
 43                 }
 44                 }
 45                 $mut_num ++;
 46             }else{
 47                 print"$_\n";
 48                 $mut_other ++;
 49             }
 50     }
 51 }
 52
 53 print OUT "The Chr:Chr$chr\n";
 54 print OUT "Total number of snp and indel:$mut_num\n";
 55 print OUT "Total number of the other mutations:$mut_other\n";
 56
 57 for my$type(keys %hash){
 58     print OUT "The mutation type:$type\n"59     my $total = sum values %{$hash{$type}};
 60     for my$gt(keys %{$hash{$type}}){
 61         my $percent_gt = $hash{$type}{$gt}/$total;
 62         my $all = sprintf "%s\t%d\t%0.3f","$gt","$hash{$type}{$gt}","$percent_gt";
 63         print OUT "$all\n";
 64     }
 65 }
 66 print OUT "########################################\n";
 67 close VCF;
 68 close OUT;


苛求远离完美
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2017-1-7 15:14:32 | 显示全部楼层
这个是linux脚本:
[Bash shell] 纯文本查看 复制代码
 1 #!/bin/bash
  2
  3 perl=/share/Data01/huizhen/data/1k_seq/count_mutation.pl
  4
  5 for i in $(cd $1 | find $1 -name "*vcf.gz")
  6     do
  7         if echo $i | grep -q "chrY"
  8             then
  9             :
 10             elif echo $i | grep -q "chrMT"
 11             then
 12             :
 13             elif echo $i | grep -q "wgs"
 14             then
 15             :
 16             else
 17             gunzip -c $i > temp
 18             perl $perl temp all.out
 19             rm -rf temp
 20             fi
 21     done

苛求远离完美
回复 支持 反对

使用道具 举报

11

主题

54

帖子

242

积分

中级会员

Rank: 3Rank: 3

积分
242
QQ
发表于 2017-1-7 15:18:56 | 显示全部楼层
为啥autosome和X突变类型有:SNP,INDEL,SV(包含各种);Y染色体是MNP,SNP,INDEL和CNV;MT包含SNP,IDNEL和MNP!!!这里的计算没有Y和MT以及wgs!
苛求远离完美
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-16 09:13 , Processed in 0.043530 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.