搜索
查看: 4153|回复: 4

跟据CDS突变位置快速查找rs号

[复制链接]

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
发表于 2016-12-18 21:50:43 | 显示全部楼层 |阅读模式
    最近有做核酸飞行时间质谱的同事想知道snp的rs号来根据软件设计引物,由于NCBI网页查找很繁琐,所以看我这边有没有快速查找的方法。但是很多时候,我们没有rsID号,只有CDS区域的碱基突变信息。所以我们的重点是跟据CDS突变位置快速查找rs号!
    我根据我以往NCBI查找rs号的经验,根据所在染色体以及染色体的位置来查找的;例如已知乙醛脱氢酶ALDH2基因的一个突变位于染色体12,位置111803962(GRCh38.p7),那么我就会在NCBI的snp搜索栏输入12:111803962,回车后就会查到该rs号为rs671 (这是我常用的搜索方式,如果大家有更好的可以告诉我,谢谢!)。所以我的想法是根据突变信息(例如:c.233T>A;且常用的就是定位到基因的CDS序列位置),来查找它在染色体上的位置,然后根据染色体的位置搜索dbsnp数据库,来查找rs号!

    好了,有了思路,就要知行合一了

    首先下载GRCh38.P7的gtf文件,在enseml网站下载目前NCBI使用的参考序列对应的基因注释版本:ftp://ftp.ensembl.org/pub/releas ... ns.GRCh38.86.gtf.gz(#!genome-build GRCh38.p7,同时gtf格式可以在该下载目录下的README文件中了解到),对应gtf信息。这样我们就可以在该文件中查找一个基因的染色体位置信息了。

    这里面针对一个基因有很多转录本,根据我以往的经验,往往我只会查找001号的转录本(例如:ALDH2-001),因为该版本往往是注释信息和链接信息最多的版本,当然使用也是最多的了!

    然后我会在该转录版本下搜索出该出该基因的CDS信息和各序列起始位置!例如ALDH2-001,转录本ID为:ENST00000261733,对应的CDS名称为:CCDS9155,这样我就根据CDS所在的外显子顺序进行CDS的先后排序,然后在做出CDS碱基位置和所在染色体位置的hash关系,那么我有了突变信息就可以跟据这个关系查到所在的染色体关系了!需要注意的是我们搜索基因所处的正负链情况要区别对待,不管怎样主要知道转录时5端到3端就行了,不管正负链情况,那么我们再排列负链基因时,CDS序列时就要根据gtf文件的排序作reverse处理了(因为单个CDS的起始位置是从小到大的,那么负链基因的转录顺序就是从大到小的位置顺序来完成的,所以此处要颠倒单个CDS中的起始顺序,但是处在不同外显子中的顺序是不需要的!)。如果比对下NCBI的CCDS序列时,会发现,我们再gtf搜索到的CDS起始位置会比NCBI的CCDS长度少3个碱基,那是因为NCBI没有吧终止密码算上去的缘故(例如:ALDH2基因的gtf起始位置为:111766983-111809572,而NCBI的CCDS起始位置为:[size=13.3333px]111766983-[size=13.3333px]111809575,这里的终止密码位置在gtf文件中也可以查到:[size=13.3333px]111809573-111809575)!此时我就可以根据突变信息查到突变所在的染色体位置了,下一步是去dbsnp搜索该位置的rs号了!

    dbsnp的bed文件下载地址(我所选的)为:ftp://ftp.ncbi.nlm.nih.gov/snp/o ... _b149_GRCh38p7/BED/,这里是根据染色体分开来记录的,我使用了wget写了脚本来下载了所有的文件,然后合并到一起!关于bed格式说明见该地址下文件:ftp://ftp.ncbi.nlm.nih.gov/snp/s ... _Mapping.README.txt。这样我就可以写一个简单脚本去bed文件中查找rs号了,对了我选择的是突变为的end位置,具体为啥不选start位置,我也没搞清楚,求路过的小盆友给点帮助
chrom        chr1         Reference sequence chromosome or scaffold; Exception: (see NOTE #5)        ExchangeSet/Rs/Assembly/Component/@chromosome                        
  chromStart        10144        Start position in chrom        ExchangeSet/Rs/Assembly/Component/@start                          chromEnd        10145        End position in chrom        ExchangeSet/Rs/Assembly/Component/@end                             name        rs144773400        dbSNP Reference SNP (rs) identifier        ExchangeSet/Rs/@rsId)!

    另外我还下载了GRCh38.P7的refseq,这样我还查到了该突变前后碱基情况,在输出前后突变信息的时候,也要考虑基因所在的正负链问题了!
    打了折磨多字,我总计下我的思路:
1,根GTF文件查到所在染色体位置;
2,在dbsnp文件中,根据染色体所在位置,查到rs号;
3,根据参考基因组输出突变碱基前后的序列情况(当然这里是基因组的序列)!



   最后我贴上我的代码,大家尽管提意见,我是菜鸟选手,才学生信5个月,还请多多指正!





上一篇:不能成功安装Rtools
下一篇:TCGA中提取lncRNA
苛求远离完美
回复

使用道具 举报

0

主题

2

帖子

24

积分

新手上路

Rank: 1

积分
24
发表于 2017-1-5 17:59:41 | 显示全部楼层
我对编程的事情不是很清楚,但其实我有个简单的办法:
http://bioinformatics.mdanderson.org/transvarweb/
TransVar可以根据变异在cDNA, gDNA, 蛋白质上的位置来提供注释:
选择cDNA作为提交的变异所在序列,选择hg19,选择用Ensembl来注释,然后按要求提供变异的标识:PIK3CA:c.1633G>A
最后就能在结果里看到需要的信息了。
回复 支持 1 反对 0

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
 楼主| 发表于 2016-12-18 21:53:11 | 显示全部楼层
[Perl] 纯文本查看 复制代码
#!/home/huizhen/bin/perl/bin

use strict;
use warnings;
#use Smart::Comments;
#use List::Util qw(sum max min);

die"perl $0 <gene> <snp_ccds_position>:$!" unless @ARGV == 2;
my $gene = shift @ARGV;
my $pos = shift @ARGV;

open GTF,"<","/share/Data01/huizhen/reference/NCBI/Homo_sapiens.GRCh38.86.gtf" or die "cant open:$!";
if(-e "$gene.info.txt"){
		die "The output file exists:$!";
}else{
		open OUT,">","$gene.info.txt";
}
my %cds;
my $strand;
my ($chr_of_snp,$exon_of_snp,$pos_of_snp);
my (@gene,@trans,@exon,@cds,@start_end,@utr);
push @gene,"The info of gene\n";
push @trans,"The info of transcript\n";
push @exon,"The info of exon\n";
push @cds,"The info of CDS\n";
push @start_end,"The info of Start_end code\n";
push @utr,"The info of UTR\n";
my @temp;
while(<GTF>){
		if(/^#/){
				print OUT $_;
		}else{
				chomp;
				if(/\tgene\t.*gene_name "$gene"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @gene,"The gnen:$gene\tThe gene_id:$1\n" if ($_ =~ m/gene_id "(\w+)";/);
						push @gene,"The seqname:$info[0]\tThe source:$info[1]\tThe feature:$info[2]\n";
						push @gene,"The start:$info[3]\tThe end:$info[4]\tThe strand:$info[6]\tThe length:$len\n";
						$chr_of_snp = $info[0];
						$strand = $info[6];
				}elsif(/\ttranscript\t.*transcript_name "$gene-001";/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @trans,"The transcript_id:$1\tThe transcript_name:$2\n" if($_ =~ m/transcript_id "(\w+)";.*transcript_name "(\w+\-?\w+)";/);
						push @trans,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe length:$len\n";
				}elsif(/\texon\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @exon,"The exon_number:$1\tThe exon_id:$2\n" if($_=~ m/exon_number "(\w+)";.*exon_id "(\w+)";/);
						push @exon,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe length:$len\n";
				}elsif(/\tCDS\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @cds,"The ccds_id:$2\tThe exon_number:$1\n" if($_ =~ m/exon_number "(\w+)";.*ccds_id "(\w+)";/);
						push @cds,"The feateur:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe frame:$info[7]\tThe length:$len\n";
						$cds{$1} = [] if not exists $cds{$1};
						######$info[3] += $info[7]; ===>>> NCBI web include the frame and start code , but not the stop code!
						######==========================>> The gtf file include the start and stop codes !
						push @{$cds{$1}},($info[3]..$info[4]);
						push @temp,$len;
				}elsif(/\tstart_codon\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @start_end,"The transcript_id:$1\tThe exon_number:$2\n" if($_ =~ m/transcript_id "(\w+)";.*exon_number "(\w+)";/);
						push @start_end,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe frame:$info[7]\tThe length:$len\n";
				}elsif(/\tstop_codon\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @start_end,"The exon-number:$1\n" if($_ =~ m/exon_number "(\w+)";/);
						push @start_end,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe frame:$info[7]\tThe length:$len\n";
				}elsif(/\tfive_prime_utr\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @utr,"The gene_id:$1\tThe transcript_id:$2\n" if($_ =~ m/gene_id "(\w+)";.*transcript_id "(\w+)";/);
						push @utr,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe length:$len\n";
				}elsif(/\tthree_prime_utr\t.*transcript_name "$gene-001"/){
						my @info = split;
						my $len = $info[4] - $info[3] + 1;
						push @utr,"The feature:$info[2]\tThe start:$info[3]\tThe end:$info[4]\tThe length:$len\n";
				}
		}
}
close GTF;
my @total = (@gene,@trans,@start_end,@exon,@cds,@utr);
for(@total){
		print OUT $_;
}
close OUT;

my %seq;
my ($position_f,$position_r)=(0,0);
for my $num(sort {$a<=>$b}keys %cds){
		if($strand eq "+"){
				for(@{$cds{$num}}){
				    $seq{$position_f} = $_;
				    $position_f ++;
				}
		}else{
			  for(reverse @{$cds{$num}}){
						$seq{$position_r} = $_;
						$position_r ++;
				}
		}
}

$pos_of_snp = $seq{$pos};
print"The snp location of Genome:chr$chr_of_snp\tThe strand:$strand\n";

for my $num(sort {$a<=>$b} keys %cds){
		for(@{$cds{$num}}){
				$exon_of_snp = $num if ($pos_of_snp == $_);
		}
}
print"The snp locaton of Exon:$exon_of_snp\n";
print"The snp position of Chromosome:$pos_of_snp\n";

open BED,"<","/share/Data01/huizhen/reference/NCBI/bed_chr_all.bed" or die "cant open:$!";
my $rs ||= "None";
while(<BED>){
		if(/$chr_of_snp\t\d+\t\b$pos_of_snp\b/i){
		    $rs = (split /\t/,$_)[3];
		    last;
		}
}
close BED;
print"The snp of rs:$rs\n";
my $genome_chr = join "",("chr_","$chr_of_snp");
open SEQ,"<","/share/Data01/huizhen/reference/NCBI/GRCh38.p7/$genome_chr" or die"cant open:$!";
my $seq_of_chr;
while(<SEQ>){
		chomp;
		$seq_of_chr .= $_ unless ($_ =~ m/^>/);
}
my $seq_before = substr($seq_of_chr,$pos_of_snp - 26,25);
my $seq_mut = substr($seq_of_chr,$pos_of_snp,1);
my $seq_after = substr($seq_of_chr,$pos_of_snp,25);
print"The bases before and after the mutation:\n";
if($strand eq "+"){
    print"$seq_before     $seq_mut     $seq_after\n";
}else{
		$seq_before =~ tr/atcgATCG/tagcTAGC/;
		$seq_before = reverse $seq_before;
		$seq_mut =~ tr/atcgATCG/tagcTAGC/;
		$seq_mut = reverse $seq_mut;
		$seq_after =~ tr/atcgATCG/tagcTAGC/;
		$seq_after = reverse $seq_after;
		print"$seq_after     $seq_mut     $seq_before\n";
}
苛求远离完美
回复 支持 反对

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
 楼主| 发表于 2016-12-18 21:55:07 | 显示全部楼层
好吧,大概思路就是这样,我写的代码,还打印了该基因的信息到输出文件中,屏幕输出为rs信息和碱基序列,同时,我还对refseq根据染色体进行了拆分,便于输出突变前后的碱基序列!
苛求远离完美
回复 支持 反对

使用道具 举报

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
 楼主| 发表于 2017-1-6 13:25:22 | 显示全部楼层
jiongz 发表于 2017-1-5 17:59
我对编程的事情不是很清楚,但其实我有个简单的办法:
http://bioinformatics.mdanderson.org/transvarweb/ ...

哇塞,我看看,谢谢!
苛求远离完美
回复 支持 反对

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-11-21 10:02 , Processed in 0.034207 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.