Researching SNPs TSD Journal

Back to results.

HGMD

Sequence retrieval and comparison

Retrieve the protein sequence from HGMD by going to the sequence view and applying switch view, as described in the task description. Download the page manually and save to a file (cdna_new.php). Then execute the following command: <source lang="bash">

grep -P --only-matching "\w{3}
\w{3}" ./cdna_new.php | cut -c99-101 | sed -n '$!p' | tr -d '\n' | tr '[:lower:]' '[:upper:]' |

./ThreeToOneLetterCode.py | tr -d '\n' > hgmdSequence_aa </source> Where 'ThreeToOneLetterCode.py' is the following script <source lang="python">

1. !/usr/bin/env python

from Bio.PDB import to_one_letter_code as one_letter import sys

seq = sys.stdin.readline()

for aaa in range(0, len(seq), 3) :

 print one_letter[seq[aaa:aaa+3]]

</source> To check whether the protein sequence from HGMD is the same than the one in Uniprot perform the following operations: <source lang="bash"> curl -s http://www.uniprot.org/uniprot/P06865.fasta | sed '1d' | tr -d '\n' > HEXA_HUMAN #Get Uniprot sequence without header perl -p -e 's|(.)|$1\n|g' HEXA_HUMAN > temp1 #Insert linebreak after every character in both sequences, for easy diff'ing perl -p -e 's|(.)|$1\n|g' hgmdSequence_aa > temp2 diff temp1 temp2 </source> This yields that position 436 differs with an Ile in Uniprot and a Val in HGMD.

Mutation retrieval and parsing

Navigate to the subsite with the missense/nonsense mutations and manually save it to 'hgmd_hexa_misssense_nonsense.php'. Then apply the following command: <source lang="bash"> parseHGMDMutations.pl hgmd_hexa_misssense_nonsense.php > missenseSNPHGMD_NONPROF.txt perl parseHGMDMutations.pl -h ./hgmd_hexa_misssense_nonsense_FULL.php -p > missenseSNPHGMD.txt #If given access to the professional version

1. Afterwards manually remove non-unique entries, if present

</source> Where 'parseHGMDMutations.pl' is the following script: <source lang="perl">

1. !/usr/bin/perl

use strict; use warnings; use sigtrap; use autodie; use Bio::SeqUtils;

my $professional = 0; my $hgmd_page;

use Getopt::Long;

   my $data = "file.dat";
   my $length = 24;
   my $verbose;
   my $result = GetOptions (
   "h=s" => \$hgmd_page, # string
   "p" => \$professional); # flag

if(!defined $hgmd_page) {

   die "Define path to hgmd PHP file\n";

}

my $fh; open($fh, '<', $hgmd_page); while(my $line = <$fh>) {

   if($professional)
   {

# cATG-GTGMet1Valc.1A>Gp.M1VTay-Sachs disease<a href='http://www.ncbi.nlm.nih.gov/sites/entrez?cmd=Retrieve&db=PubMed&list_uids=1532289&dopt=Abstract' target=ref>Mules (1992) Am J Hum Genet 50, 834</a> <img src='bar/dm.gif' title='Disease causing mutation' alt='Disease causing mutation'> <a href='http://www.ncbi.nlm.nih.gov/sites/entrez?db=snp&cmd=search&term=rs121907965' target='dbSNP-entry'><img src='bar/dbsnp.gif' border='0' title='dbSNP ID: rs121907965' alt='dbSNP ID: rs121907965'></a> <img src='bar/genomic.gif' title='Coordinate Chr 15: 72668313' alt='Coordinate Chr 15: 72668313'> if($line =~ m|\w+-\w+(\w{3})(\d+)(\w{3}).+p\..+(.+?).+|)

       {
           my $id = ;
           my $aa_wt = $1;
           my $aa_mt = $3;
           my $pos = $2;
           my $pheno = $4;
           
           print $Bio::SeqUtils::ONECODE{$aa_wt} . $pos . $Bio::SeqUtils::ONECODE{$aa_mt} . "\tHGMD\t$id\t$pheno\tmissense\n";
       }
   }
   else
   {
       #This pattern matching will not find nonsense mutations, therefore already implicitly excluding these

if($line =~ m|<a name="\w{8}">(\w{8})</a>\w+-\w+(\w{3})-(\w{3})(\d+)<a href="http://www.biobase-international.com/product/hgmd"><img src=".+" alt="BioBase" border="0"></a>(.+).+|)

       {
           my $id = $1;
           my $aa_wt = $2;
           my $aa_mt = $3;
           my $pos = $4;
           my $pheno = $5;
           
           print $Bio::SeqUtils::ONECODE{$aa_wt} . $pos . $Bio::SeqUtils::ONECODE{$aa_mt} . "\tHGMD\t$id\t$pheno\tmissense\n";
       }
   }

} close($fh); </source> To find the intersection and unique elements between the free and professional version of HGMD use the general mechanism exemplified for SNPedia.

dbSNP

At first the SNPs were queried with "synonymous-codon"[Function_Class] AND HEXA[GENE] AND "human"[ORGN] AND "snp"[SNP_CLASS] in dbSNP. Those results were saved and then the mutations extracted with

bash getdbSNP.sh entrez.html

It is important to note that in the dbSNP results the rs ids can be redundant. getdbSNP.sh: <source lang="bash">

1. !/bin/bash

${1?"Usage: $0 textfile"}

website=$1

rm -f Mut_$website

for rs in `egrep -o ">rs[0-9]+<" $website |uniq | egrep -o "rs[0-9]+"`; do echo $rs >> Mut_$website curl -s "http://www.ncbi.nlm.nih.gov/projects/SNP/snp_gene.cgi?rs=$rs" | egrep "aaCode|proteinPos" | cut -d \" -f 4 >> Mut_$website echo "---" >> Mut_$website done </source>

The mutations were then checked with the reference sequence (!) and brought into the right format:

perl formatSNPs.pl -i Mut_entrez.html -f P06865.fasta

<source lang="perl">

1. !usr/bin/perl

use Getopt::Long; our $opt_f = (); our $opt_d = "dbSNP"; sub usage {

       printf "Usage: $0 [flags] file...\n";
       printf "flags:\n";

       printf "    -i Infile   raw mutations from getdbSNP.sh \n";
       printf "    -f  Referenzfile fastafile containing referenz sequenz ";
       printf "    -missense  Missesne mutations are considered, default = synonymous ";
       printf "    -d name of database to write in output file, default = dbSNP ";
       exit 1;

} GetOptions( "i=s", "f=s", "missense", "d=s" ); usage() if !( $opt_i && $opt_f );

open( REF, "<$opt_f" ) or die "Referenzfile not found";

@reference = (); while (<REF>){ chomp; @line = split( , $_ ); if ($line[0] ne ">") { @reference = (@reference,@line); } }

open (IN, "<$opt_i") or die "Inputfile not found"; %muts = (); my $rs; while(<IN>){ chomp; if ($_ eq "---"){ if (!$opt_missense){ if ($reference[$pos-1] eq $aa){ print "$aa$pos$aa\t$opt_d\t$rs\t\t". (($opt_missense) ? 'missense' : 'synonymous') ."\n"; } else { print "! not matching reference $reference[$pos-1]: $aa$pos$aa\n"; } } else { foreach my $key ( keys %muts ){ if ($key ne $aa){ if ($reference[$pos-1] eq $aa){ print "$aa$pos$key\t$opt_d\t$rs\t\t". (($opt_missense) ? 'missense' : 'synonymous') . "\n"; } else { print "! not matching reference $reference[$pos-1]: $aa$pos$key\n"; } } } } $count=0; $aa=""; %muts = ();

       $rs = ;

} elsif ($count==0){ $rs=$_; $count++; } elsif ($count == 1)

   {
     $pos = $_;
     $count++;
   }

elsif ($count>=2){ $muts{$_}=1; $aa=$_; $count++; } } </source>

The same procedure can be performed with the missense Mutations from dbSNP with the dbSNP search "missense"[Function_Class] AND HEXA[GENE] AND "human"[ORGN] AND "snp"[SNP_CLASS] followed by

bash getdbSNP.sh entrezMissense.html
perl formatSNPs.pl -i Mut_entrezMissense.html -f P06865.fasta -missense

SNPdbe

The results are provided on the webpage in a tabular format, where a lot of information can be obtained very easily just by sorting the columns (e.g. conservation scores, exp. evidence). A text file with the results snps.txt can also be downloaded directly from the webpage. Within the text file the fields are separated by "|". In the first field the mutations is displayed. To format the SNPs according to the previously used conventions SNP ORIGIN ID this script can be used.

<source lang="perl">

1. !usr/bin/perl

open(IN,"<$ARGV[0]") or die "no input file"; while (<IN>){

       $sep="";
       $iddbsnp="";
       $idswissvar="";
       $Dsep="";
       @line = split('\|', $_);
       print $line[0]."\tsnpdbe: ";
       if ($line[1] ne "N/A"){
               print "dbSNP";
               $sep=",";
               $Dsep=",";
               $iddbsnp = $line[1];
       }
       if ($line[2] ne "N/A"){
               print $sep."SwissProt";
               $sep=",";
       }
       if ($line[3] ne "N/A"){
               print $sep."SwissVar";
               $sep=",";
               $idswissvar = $line[3];
       } else { $Dsep="";}
       if ($line[4] ne "N/A"){
               print $sep."PMD";
               $sep=",";
       }
       if ($line[5] ne "N/A"){
               print $sep."1000genomes";
       }
       print "\t $iddbsnp$Dsep$idswissvar \n";

} </source>

The conservation scores are stored in the field 21 and 24, perc_wt and perc_mt:

cut -d "|" -f 1,22,25 snps.txt | sed s/\|/" "/g

OMIM

Entries with rs-identifiers

Retrieve all entries that have rs-Identifiers: <source lang="bash"> curl -s http://omim.org/allelicVariant/606869 | grep -P --only-matching "target=\"_blank\">rs\d+</a>" | sed 's/\s//' | grep -P --only-matching "rs\d+" | sort | uniq > rs_OMIM.txt </source> Get more information on the entries from dbSNP, using the same script as used for SNPedia: <source lang="bash"> ../snpedia/getSNPData.sh rs_OMIM.txt </source> and use the same script as used for dbSNP to format the final files with missense and synonymous SNPs: <source lang="bash"> perl ../dbSNP/formatSNPs.pl -i ./enriched_mis_rs_OMIM.txt -f ../dbSNP/P06865.fasta -missense -d OMIM > missenseSNPOMIM.txt perl ../dbSNP/formatSNPs.pl -i ./enriched_syn_rs_OMIM.txt -f ../dbSNP/P06865.fasta -d OMIM > synonymousSNPOMIM.txt </source> Intersection to dbSNP can be found analogous to the description for SNPedia. Note that in this case only 29 rs-identifiers are supposed to be in the intersection. However since one identifer appears two times in OMIM (it describes two mutations), there is actually no entry unqiue to OMIM. This can also be verified with the following commandline (using the same naming scheme as for SNPedia): <source lang="bash"> grep -Fxv -f dbsnpc tempc #Does not return anything </source>

Finally add the one entry that does not have an rs-Identifier but describes a valid missense variant (C58Y).

SNPedia

Retrieve all rs-identifiers from SNPedia (and already exclude one entry without protein information), using the following commandline: <source lang="bash"> curl -s http://www.snpedia.com/index.php\?title\=Special:Search\&ns0\=1\&redirs\=0\&search\=Gene%3DHEXA\&limit\=500\&offset\=0 | grep -P --only-matching "title=\"Rs\d+\"" | sed 's|title=||' | tr -d '"' | tr '[:upper:]' '[:lower:]' | grep -v 'rs121907980' > rs_SNPedia.txt </source> Then get more information on the rs-identifiers from dbSNP with the following comamndline <source lang="bash"> ./getSNPData.sh rs_SNPedia.txt </source> where getSNPData.sh ist the following script; <source lang="bash">

1. !/usr/bin/bash

in=$1 rm -f enriched_syn_$in rm -f enriched_mis_$in

for rs in `cat $in`; do

 curl -s "http://www.ncbi.nlm.nih.gov/projects/SNP/snp_gene.cgi?rs=$rs"  > temp
 if [ `egrep 'missense' ./temp |wc -l` -gt 0 ]
 then
   echo $rs >> enriched_mis_$in
   cat temp | egrep "aaCode|proteinPos" | cut -d \" -f 4 > temp2
   #In theory, do special parsing of entries that consist of two or more mutations. Luckily for our case there are none ;)
   cat temp2 >> enriched_mis_$in
   echo "---" >> enriched_mis_$in
 else
     if [ `egrep 'synon' ./temp |wc -l` -gt 0 ]
     then
       echo $rs >> enriched_syn_$in
       cat temp | egrep "aaCode|proteinPos" | cut -d \" -f 4 > temp2
       #In theory, do special parsing of entries that consist of two or more mutations. Luckily for our case there are none ;)
       cat temp2 >> enriched_syn_$in
       echo "---" >> enriched_syn_$in
     else
       echo "Ignoring $rs because neither missense nor synon"
     fi
   fi

done

rm -f ./temp rm -f ./temp2 </source> With the output of this, we can use the script TODO from above to get the final formatted SNPs: <source lang="bash"> perl ../dbSNP/formatSNPs.pl -i ./enriched_mis_rs_SNPedia.txt -f ../dbSNP/P06865.fasta -missense -d SNPedia > missenseSNPSNPedia.txt perl ../dbSNP/formatSNPs.pl -i ./enriched_syn_rs_SNPedia.txt -f ../dbSNP/P06865.fasta -d SNPedia > > synonymousSNPSNPedia.txt </source>

Overlap to dbSNP

<source lang="bash"> cut -f 3 missenseSNPSNPedia.txt > temp1 cut -f 3 synonymousSNPSNPedia.txt > temp2 cat temp1 temp2 > tempc cat ../dbSNP/synonymous_rsids.txt ../dbSNP/missense_rsids.txt > dbsnpc grep -Fx -f tempc dbsnpc | wc -l #Yields: 30 </source>

Mutation map

Create a complete map of mutations by combining all previous files and then call the master script that creates all source files for plotting: <source lang="bash"> cat ../dbSNP/missenseSNPdbSNP.txt ../dbSNP/synonymousSNPdbSNP.txt ../OMIM/missenseSNPOMIM.txt ../OMIM/synonymousSNPOMIM.txt ../HGMD/missenseSNPHGMD.txt ../snpdbe/mutSnpdbe.txt ../snpedia/missenseSNPSNPedia.txt ../snpedia/synonymousSNPSNPedia.txt > combinedSNP.txt perl createMap.pl combinedSNP.txt > mutationMap #More files are created in the background </source> createMap.pl

Plotting

Wildtype and Mutant Amino Acids

<source lang="bash"> library(gplots) #Needs caTools, bitops and a few more. Gives us: barplot2

clwdaxis <- 20 ccaxis <- 2 ccmain <- 1.7 cbarnames <- 2 cclab <- 2.7

d <- read.table("aaPropensities")

colors <- c("DarkGreen", "red")

png(width=1000, height=900, "TSD SNP aaPropensities.png") par(mar=par()$mar+c(0,1,-1,0)) barplot2(as.matrix(d), beside=TRUE, col=colors, cex.axis=ccaxis, cex.main=3.0, cex.lab=cclab, cex.names=cbarnames, main="Frequency of wildtype and mutant amino acids", ylab="# of occurences",ylim=c(0,20)) legend("topleft", c("Wildtype", "Mutant"),fill=colors, cex=2.5, inset=0.05) dev.off() </source>

Databases per positon

TODO

Pymol

Based on the script by Fabry-Disease. Execute, save output to script.pml and execute with python -cq script.pml <source lang="bash">

1. !/bin/bash

cat <<EOT reinitialize load referenz_close_lessshallow.pse, referenz bg_color white set ray_opaque_background, off

1. 1. 1. non disease causing###

unset resi for resi in `cat resi_nsynndisonly`;do echo "select non-disease, resi $resi, merge = 1" done if -n $resi ; then cat <<EOT color blue, non-disease disable non-disease

EOT fi

1. 1. 1. disease###

unset resi for resi in `cat resi_nsyndisonly`;do echo "select disease, resi $resi, merge = 1" done if -n $resi ; then cat <<EOT color red, disease disable disease EOT fi

EOT

1. 1. 1. silent###

unset resi for resi in `cat resi_synonly`;do

   echo "select silent, resi $resi, merge = 1"

done

if -n $resi ; then

   cat <<EOT

color slate, silent disable non-disease

set ray_trace_mode, 0; # color raytrace, no outlines set antialias, 2; #High anti-aliasing

1. How do I crank-up the glossyness of rendered atoms?

set spec_power = 200; # in range 20 - 250 set spec_refl=1.5; # in range 1.0 - 2.0

set orthoscopic, off; #perspective set ray_shadows,0; # turn off shadows ray 1450,900; save highStructure.png EOT fi </source>