tag:blogger.com,1999:blog-27357508530961824912024-03-04T23:14:44.422-08:00genetics for funUnknownnoreply@blogger.comBlogger15125tag:blogger.com,1999:blog-2735750853096182491.post-67090943356455202832016-12-18T21:16:00.003-08:002017-01-30T20:02:00.505-08:00What can we learn from one Y chromosome?Some time ago I became interested in chromosome Y. It all started from my maternal grandfather. He was born in Lauria, a small town in the south of Italy. After the civil records for Lauria became available online, it was pretty easy to reconstruct his family tree all the way to the 18th century. However, as shown in the following diagram, one small piece was missing (at the top right):
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEinKNQK0SbhvLshXZ0wuA_InKIgMJkCIa-_06dJLKoWq3Aa_UqWzDM-9_YazkXp_t1wtu9ugKqeF1zGOHUknP7awxiSLA8XLvERLjU8Zee2zZMgNLKcN95Z_eF_emsFwmoG-NhRkoNcYdg/s1600/grandpa.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEinKNQK0SbhvLshXZ0wuA_InKIgMJkCIa-_06dJLKoWq3Aa_UqWzDM-9_YazkXp_t1wtu9ugKqeF1zGOHUknP7awxiSLA8XLvERLjU8Zee2zZMgNLKcN95Z_eF_emsFwmoG-NhRkoNcYdg/s1600/grandpa.png" width="560" /> </a>
<br />
My grandfather's great-grandfather on the direct paternal line was listed as a child of unknown father. His mother got <a href="http://antenati.san.beniculturali.it/v/Archivio+di+Stato+di+Potenza/Stato+civile+napoleonico/Lauria/Matrimoni/1812/007694080_01373.jpg.html">married in 1812</a> to someone who subsequently disappeared from the record, and <a href="http://antenati.san.beniculturali.it/v/Archivio+di+Stato+di+Potenza/Stato+civile+della+restaurazione/Lauria/Nati/1832/007693350_01063.jpg.html">gave birth in 1832</a> to a child as a single mother. While at first this was puzzling, and unusual for the times, it did seem to fit a particular historical context. As a child she had lost her father in 1806 during the <a href="https://it.wikipedia.org/wiki/Massacro_di_Lauria">massacre of Lauria</a> that saw the population of the town slaughtered by the French soldiers under general André Masséna, as a punishment for having supported the Bourbon kings during the Napoleon invasion. The French presence brought an increase of secular behaviours among which a significant increase of single mothers (<a href="https://books.google.com/books?id=3jWq7N9_cmMC&pg=PA105">ref</a>).
<br />
<br />
Chances are that my 4th great-grandfather on the direct paternal line of my grandfather was another local person from Lauria, who decided to never recognize his own biological child who died in 1911 still with the father reported as unknown. Whatever happened, it seemed one of those maybe irrelevant mysteries doomed to be forgotten in time. But then I thought that there was a chance, albeit slight, to shed some light over what happened. Whoever this unknown father who conceived my 3rd great grandfather in 1831 was, he must have carried the same Y chromosome my maternal grandfather carries.
<br />
<br />
Now, my maternal grandfather does not have sons, does not have brothers, and all of his paternal uncles either left no sons or they moved to Colombia in the 20th century, eventually losing touch with my family. He truly is the only person in my family that carries the Y chromosome that could unlock the mystery. Last year I struck a deal with him promising that, in exchange of some of his DNA, one day I will figure out where his patrilineal line came from and therefore what his real family name should have been if it indeed followed the patrilineal line. How does this work?
<br />
<br />
While AncestryDNA and 23andMe offer the most comprehensive autosomal tests, there is only one company, FamilyTreeDNA, which performs in depth Y chromosome analyses. The Y chromosome is a special molecule. It never recombines and it gets passed father to son as a single block. Excluding the tips, it is truly the largest molecule to withstand the test of time almost unscathed. But every Y chromosome is unique. In fact, it is so large, that we expect to find, more often than not, a single point mutation difference between the Y chromosome of a son and the Y chromosome of a father in the 21.3 megabases of DNA sequence of the chromosome Y that we can currently assay (see <a href="http://dx.doi.org/10.1038/ng.3171">here</a>). And that holds only for point mutations.
<br />
<br />
Before next-generation sequencing became as mainstream as it is now, most analyses of the Y chromosome relied on measurements of Y chromosome short tandem repeats (STRs). STRs are small segments of repetitive DNA that, due to their repetitive nature, are prone to slippage when copied and this results in their length to change. Although we can almost say that STR analyses are falling out of fashion a bit, in the last 20 years tens of thousands of Y chromosomes have had some of their STRs analyzed. There are 100s of STRs on the Y chromosome but those that have been analyzed before the advent of next-generation sequencing have the most value when you want to compare to available databases. Currently, the cheapest way to get the largest panel is to go for the Y111 DNA test with FamilyTreeDNA. This panel of STRs is such that you can observe almost with 30% chance at least one of them to change length in one father-son transmission. This means that Y chromosomes of direct relatives might actually look alike with a Y111 test, but other than that it is almost guaranteed that every Y chromosome will get its own unique fingerprint with the test and the number of differences might give an idea, albeit vague, of how distantly related two Y chromosomes are.
<br />
<br />
So I got my grandpa's Y chromosome tested. The results are just a <a href="http://www.ysearch.org/search_view.asp?viewuid=6FD3N">list of numbers</a>, measuring the length of each STR tested, but they can tell a lot. The <a href="http://www.ysearch.org/search_view.asp?viewuid=DN4GP">closest match</a> to my grandpa's Y chromosome matches 60 out of 67 markers, which corresponds to a shared ancestor within the last 20 generations with 90% probability. This is not very helpful, but it does make sense. The match originates from Mercato San Severino, which is less than 100 miles from Lauria:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj6JVy7WusMEX7pTT9ljWG_i-7EFcfUjX8U4Sc2La7AGKm3MFQFeFs5YIkdKrhesWJKVKDUmzeWoDcIUkfBsqxZjFc1kemVvaqylnhIoJWnPc72CiOmWDlJe5t6wCvz7H4rbEE7HncMEJk/s1600/loria.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj6JVy7WusMEX7pTT9ljWG_i-7EFcfUjX8U4Sc2La7AGKm3MFQFeFs5YIkdKrhesWJKVKDUmzeWoDcIUkfBsqxZjFc1kemVvaqylnhIoJWnPc72CiOmWDlJe5t6wCvz7H4rbEE7HncMEJk/s1600/loria.png" width="560" /></a>
<br />
Furthermore the match's family name is Loria, a variant spelling for Lauria. Less interestingly, a <a href="http://www.ysearch.org/search_view.asp?viewuid=BG849">second match</a>, which matches at 54 out of 67 markers, originated from Mezzojuso in Sicily. So what have I learned? I didn't get any clue about the true family name my grandpa should have received, but I have a pretty good guess that his Y chromosome has been around the area for at least 500 years, and most likely my 4th great grandfather was someone local rather than a distant stranger passing by in town. ;-)
<br />
<br />
Although ultimately I did not find the person that passed the Y chromosome to my 3rd great grandfather, this match tells me that my grandfather's results check out and they were not the result of some weird sample swap on FamilyTreeDNA's side. Maybe one day, within my lifetime, a better closer match will show up and maybe the mystery will be solved. There is a little bit of magic involved in thinking that the answer to this riddle might be out there waiting to come out and maybe one day it will make for a unique story.
<br />
<br />
I was so impressed from what could be inferred that I decided to take the Y111 test myself. Of course, my maternal grandfather and I don't share the same Y chromosome, so in this case I was investigating a completely different paternal line. Furthermore, I am a Genovese, and it would be nice to get proof that I am not actually related to every <a href="https://en.wikipedia.org/wiki/Genovese_crime_family">Genovese</a> out there. When my <a href="http://www.ysearch.org/search_view.asp?viewuid=7KEC3">results</a> came in, disappointingly, I didn't get any good matches as my grandfather did. The best matches though did make some sense. My paternal grandfather was born in Valguarnera Caropepe in Sicily, and I was able to trace his patrilineal line back to an individual named Francesco Genovese born in Piazza Armerina from Tommaso Genovese and who got <a href="http://familysearch.org/pal:/MM9.3.1/TH-266-11119-159465-60">married in 1784</a> in Valguarnera Caropepe. The four best matches I have found (<a href="http://www.ysearch.org/search_view.asp?viewuid=VEP8S">1</a>, <a href="http://www.ysearch.org/search_view.asp?viewuid=FMTK3">2</a>, <a href="http://www.ysearch.org/search_view.asp?viewuid=U6QBJ">3</a>, <a href="http://www.ysearch.org/search_view.asp?viewuid=EBTAM">4</a>) originate from Palermo, Trapani, and Caltagirone, three cities in Sicily, the last one being less than 20 miles from Piazza Armerina:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4YM5OYJTci2ZEwzKm56MNWZTjQb4k2RBXiylA909DiTZkdrh4uMV79Rb67ntqdtorQbw0WmMkqgn-URcF5wSmTDY9jrF0se4VMvvzQsLkNq5x68MhFsdZBGs4xLQRfogiUxoqNbJe-Zk/s1600/genovese.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj4YM5OYJTci2ZEwzKm56MNWZTjQb4k2RBXiylA909DiTZkdrh4uMV79Rb67ntqdtorQbw0WmMkqgn-URcF5wSmTDY9jrF0se4VMvvzQsLkNq5x68MhFsdZBGs4xLQRfogiUxoqNbJe-Zk/s1600/genovese.png" width="560" /></a>
<br />
However, these best matches still mismatch my Y chromosome at so many STRs that the time to the most recent common ancestor is difficult to pinpoint. It seems very possible that he might have lived more than 2,000 years ago. At least that tells me that my Y chromosome has been in Sicily for a very long time and it is cool that I can infer that.
<br />
<br />
What happens if we want to know about a more distant history of our Y chromosomes? When and from where did my grandfather Y chromosome and my own Y chromosome get to Italy? It turns out that there is quite some literature on the topic. Whole genome sequence analysis of 1,000 Genomes project individuals paints a picture of a sudden expansions of Y chromosome lineages during the last 10,000 years ago, together with an expansion of lineages around 50,000 years, that is, at the time when the first humans migrated out of Africa (<a href="http://dx.doi.org/10.1038/ng.3559">ref</a>):
<br />
<a href="http://www.nature.com/ng/journal/v48/n6/fig_tab/ng.3559_F2.html" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh3JE6NllFLyMTCFW1H28r2Dy-9K2-7_7KAQZtRHB91TnmwP7Q1-VqipDq0DechAc8mIjPONHLgCmmbZ46n5tfPvrD1fuv4BoA0gnKV5NYYhKYJo25BhCPmSrLOw9g5X4Z8WgEmm3A5sAo/s1600/ng.3559-F2.jpg" width="560<br" /> </a>
<br />
My Y chromosome belongs to the <a href="https://en.wikipedia.org/wiki/Haplogroup_J-M172">J2-M172</a> lineage (or more specifically J2a-L70) and my grandfather's Y chromosome belongs to the <a href="https://en.wikipedia.org/wiki/Haplogroup_R1b">R1b-M343</a> lineage (or more specifically <a href="https://en.wikipedia.org/wiki/Haplogroup_R1b#R1b1a1a2_.28R-M269.29">R1b-M269</a>), which means that they don't share an evolutionary past in the last 50,000 years.
<br />
<br />
We know a lot about the history of Y chromosomes up to the definition of each lineage. From the Y chromosome that gave origin to all Y chromosomes in people alive today, we can map many of the migrations that established the geography of nowadays geographical distribution:
<br />
<br />
<a href="https://upload.wikimedia.org/wikipedia/commons/c/ca/World_Map_of_Y-DNA_Haplogroups.png" imageanchor="1"><img border="0" src="https://upload.wikimedia.org/wikipedia/commons/c/ca/World_Map_of_Y-DNA_Haplogroups.png" width="560" /></a>
<br />
<br />
Understanding what happened within each lineage, that is, assigning a specific subclade to a Y chromosome, can be extremely challenging as the coming of agriculture brought an explosion of each lineage that survived and this resulted in some level of saturation of the STR space. Part of the reason is the massive bottleneck Y chromosomes experienced at the onset of agricolture. Around 8-4 thousands years ago the effective population size of Y chromosomes plummeted in non-African populations, possibly indicating a change in social structures that increased male variance in offspring number (<a href="http://dx.doi.org/10.1101/gr.186684.114">ref</a>):
<br />
<a href="http://genome.cshlp.org/content/25/4/459/F2.expansion.html" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjA0OuJLleTVDA-pmqDFqfVXJWI7l893PMNvkqrigStknTepSOvXATIC6ZJ5hqpCn2bqPewyn7ug0SShebPyDEGg5PG9eQbx9dX-zhR-Fm378Vlir7rmyyiibOz3W4R8MfTVyblGrlI4bo/s1600/F2.large.jpg" width="560" /></a>
<br />
Whatever the reason, the consequences are quite evident in today's Y chromosome pool. Comparing your own Y chromosome to many other Y chromosomes out there alone can make this effect quite evident. Anybody can do this. In fact, FamilyTreeDNA hosts a large number of projects attempting to cluster Y chromosomes, and these projects include thousands of Y chromosome STR haplotypes available for download, provided you have a FamilyTreeDNA account and can run a <a href="https://github.com/freeseek/getftdna">script</a> that will download all the publicly available data. What happens if I compare the number of mismatched STRs (excluding multi-copy and fast-changing STRs, leaving 86 STRs) between my Y chromosome and every other Y chromosome analyzed with the Y111 kit? Here is the figure:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaVwvEGBw3NovCt8l5rlQORCajlVfFVv-bzgfW5WwbrLY67OZekQGNa3mZY9RaYijm7Hhq1JKsdXu-2g6kN8tg5PPqL9JuHlgJeYefGjIQOAfk_QhRl-Gy8dJ01PH-oGCqf4uxOXpCRMo/s1600/333499.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiaVwvEGBw3NovCt8l5rlQORCajlVfFVv-bzgfW5WwbrLY67OZekQGNa3mZY9RaYijm7Hhq1JKsdXu-2g6kN8tg5PPqL9JuHlgJeYefGjIQOAfk_QhRl-Gy8dJ01PH-oGCqf4uxOXpCRMo/s1600/333499.png" width="560" /></a>
<br />
While FamilyTreeDNA simply assigns my Y chromosome to lineage <a href="https://yfull.com/tree/J-M172/">J2-M172</a>, which formed about 27,900 years ago, the two STR values DYS445=6 and DYS391=9 clearly assign it to <a href="https://yfull.com/tree/J-L70/">J2a-L70</a>, a subclade that formed approximately 6,900 years ago and is now spread <a href="https://www.google.com/maps/d/u/0/viewer?mid=1XMZkpib63UEQMfyOLuiGo0wkO_E">all over Europe</a>. This is a subclade of clade <a href="https://yfull.com/tree/J-L24/">J2a-L24</a>, which formed around 13,800 years ago. The presence of my Y chromosome in Sicily might be explained by the colonization of the Magna Graecia by the Greeks, but a few other scenarios are also possible. The same picture for my grandfather's Y chromosome looks substantially different:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhCdD8hD9_KkdNLyONyz-LcIO-wKkJaDIa9RQ3sa75inZOkcdjgJgl7oth8BLlgef_xckiTBxB4a9Ekuc6ymc9Ek3y0VNPQp0ihT-V1xNcVwRbW-vBjdJQY-DYdSo8mXGSrfplhwMFbIJo/s1600/456438.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEhCdD8hD9_KkdNLyONyz-LcIO-wKkJaDIa9RQ3sa75inZOkcdjgJgl7oth8BLlgef_xckiTBxB4a9Ekuc6ymc9Ek3y0VNPQp0ihT-V1xNcVwRbW-vBjdJQY-DYdSo8mXGSrfplhwMFbIJo/s1600/456438.png" width="560" /> </a>
<br />
FamilyTreeDNA assigns my granpa's Y chromosome to lineage <a href="https://yfull.com/tree/R-M269/">R1b-M269</a>, which formed about 6,500 years ago. While understanding subclades within this lineage is challenging, as the number of subclades somewhat saturates the STR space, a plausible candidate is <a href="https://yfull.com/tree/R-P312/">R1b-P312</a>, a subclade that formed approximately 4,600 years ago and is now common in Italy. In this second figure we also observe a second mode mostly corresponding to Y chromosome haplotypes from the <a href="https://yfull.com/tree/R-M417/">R1a-M417</a> cluster, a sister lineage that formed around 5,500 years ago.
<br />
<br />
Finally, in the <a href="http://apol1.blogspot.com/2016/10/1000-genomes-project-phase-3-principal.html">previous post</a> I mentioned that 1,000 Genomes project individual <a href="http://catalog.coriell.org/0/Sections/Search/Sample_Detail.aspx?Ref=HG01944">HG01944</a> showed up as an outlier. While claiming that four of his grandparents were born in Peru, his autosomal DNA projected more or less in the middle of East Asians (CHB, JPT, CHS, CDX, KHV) and Peruvians (PEL):
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgsCRyNEducJ8NMIwGOPPYKevxIJurHk0KpSBOB8pQvAUh9BbOX-OAbFReOl-zJnLjn51L8a8BsYZ1Mfb68fFxEcqb1da_LHKkHJapNC79LKHN8yv2OKEXFo_uvL2H15lg2UCKVsyhTofw/s1600/kgp.png" imageanchor="1"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgsCRyNEducJ8NMIwGOPPYKevxIJurHk0KpSBOB8pQvAUh9BbOX-OAbFReOl-zJnLjn51L8a8BsYZ1Mfb68fFxEcqb1da_LHKkHJapNC79LKHN8yv2OKEXFo_uvL2H15lg2UCKVsyhTofw/s1600/kgp.png" width="560" /> </a>
<br />
What about his Y chromosome? It belongs to <a href="https://en.wikipedia.org/wiki/Haplogroup_Q-M120">Q-M120</a> (a subclade of lineage <a href="https://en.wikipedia.org/wiki/Haplogroup_Q-M242">Q-M242</a>), a subclade separate from lineage <a href="https://en.wikipedia.org/wiki/Haplogroup_Q-M3">Q1a-M3</a> (found in more than 90% of Native American Y chromosomes). Instead, the subclade of HG01944 Y chromosome <a href="https://www.yfull.com/tree/Q-M120/">Q-M120</a> is found to some extent in East Asia, indicating that most likely the paternal grandfather of HG01944 was indeed East Asian rather than Peruvian. Another way to say it is that, while HG01944 might not be the best individual for a Peruvian reference panel, he is indeed the result of what we could call love without boundaries. ;-)Unknownnoreply@blogger.com2tag:blogger.com,1999:blog-2735750853096182491.post-24512554581308111162016-10-14T21:00:00.000-07:002016-11-03T15:36:39.244-07:001000 Genomes project phase 3 principal component analysisThe 1000 Genomes project phase 3 genotype data has been available since 2014, but I have not seen any detailed instructions for how to generate a principal component analysis plot of the 2,504 individuals for which genotype data is available. Let's fix that. First of all, you will need software to deal with the formats the data is being distributed as:
<pre>
# create your personal binary directory
mkdir -p ~/bin/
# install software needed for basic commands
sudo apt-get install git wget unzip samtools
# install latest version of bcftools
git clone --branch=develop git://github.com/samtools/bcftools.git
git clone --branch=develop git://github.com/samtools/htslib.git
cd htslib && make && cd ..
cd bcftools && make && cd ..
/bin/mv bcftools/bcftools ~/bin/
# install latest version of plink2
wget https://www.cog-genomics.org/static/bin/plink/plink_linux_x86_64_dev.zip
unzip -od ~/bin/ plink_linux_x86_64_dev.zip plink
# install latest version of GCTA
wget http://cnsgenomics.com/software/gcta/gcta_1.25.2.zip
unzip -od ~/bin/ gcta_1.25.2.zip gcta64
</pre>
Now that we have all necessary software, we need to download a copy of the human genome reference and a copy of the genotype data. These operations might take several hours, so they might be best run on a good internet connection and possibly overnight:
<pre>
# download a version of human genome reference GRCh37
wget -O- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz | gunzip > human_g1k_v37.fasta
samtools faidx human_g1k_v37.fasta
# download 1000 Genomes reference panel version 5a
mkdir -p chrs/
wget http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase3_v5a/chr{1-4,5-9,10-15,16-22X}.1kg.phase3.v5a.vcf.zip
for chrs in 1-4 5-9 10-15 16-22X; do unzip chr$chrs.1kg.phase3.v5a.vcf.zip -d chrs/; done
</pre>
The next commands are a bit technical. What we are going to is: (i) convert the downloaded VCF files into plink format using previously described <a href="http://apol1.blogspot.com/2014/11/best-practice-for-converting-vcf-files.html">best practices</a>; (ii) remove duplicate markers shamefully present in the genotype data and a few long indels that would make plink2 memory hungry; (iii) join all chromosome files into a single plink file.
<pre>
# download and convert 1000 Genomes project phase 3 reference
for chr in {1..22} X; do
tabix -f chrs/chr$chr.1kg.phase3.v5a.vcf.gz
bcftools norm -Ou -m -any chrs/chr$chr.1kg.phase3.v5a.vcf.gz |
bcftools norm -Ou -f human_g1k_v37.fasta |
bcftools annotate -Ob -x ID \
-I +'%CHROM:%POS:%REF:%ALT' |
plink --bcf /dev/stdin \
--keep-allele-order \
--vcf-idspace-to _ \
--const-fid \
--allow-extra-chr 0 \
--split-x b37 no-fail \
--make-bed \
--out chrs/kgp.chr$chr
done
# impute sex using chromosome X
plink --bfile chrs/kgp.chrX \
--keep-allele-order \
--impute-sex .95 .95 \
--make-bed \
--out chrs/kgp.chrX && \
/bin/rm chrs/kgp.chrX.{bed,bim,fam}~
# check for duplicate markers (there are 11,943 such markers, mostly on the X chromosome, unfortunately)
for chr in {{1..22},X,Y,MT}; do cut -f2 chrs/kgp.chr$chr.bim | sort | uniq -c | awk '$1>=2 {print $2}'; done > kgp.dups
# check for very long indels (there are 46 of these)
cut -f2 chrs/kgp.chr{{1..22},X,Y,MT}.bim | awk 'length($1)>=150' | sort | uniq > kgp.longindels
# generate version of each chromosome without duplicate variants
for chr in {1..22} X; do
cat kgp.{dups,longindels} |
plink --bfile chrs/kgp.chr$chr \
--keep-allele-order \
--exclude /dev/stdin \
--make-bed \
--out chrs/kgp.clean.chr$chr
done
# join all chromosomes into one
cat kgp.clean.chrX.fam > kgp.fam
cat kgp.clean.chr{{1..22},X}.bim > kgp.bim
(echo -en "\x6C\x1B\x01"; tail -qc +4 kgp.clean.chr{{1..22},X}.bed) > kgp.bed
</pre>
If you have managed to follow up to this point, the difficult part is over. We now are going to download population information and compute the principal components using plink2 and GCTA:
<pre>
# download population information
wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel
awk 'BEGIN {print "FID\tIID\tPOP"} NR>1 {print "0\t"$1"\t"$2}' integrated_call_samples_v3.20130502.ALL.panel > kgp.pop
# compute principal component analysis weights
plink --bfile kgp --maf .01 --indep 50 5 2 --out kgp
plink --bfile kgp --extract kgp.prune.in --make-grm-bin --out kgp
gcta64 --grm-bin kgp --pca 20 --out kgp --thread-num 10
(echo FID IID PC{1..20}; cat kgp.eigenvec) > kgp.pca
</pre>
The last command will generate a file which contains the first 20 principal components for the 2,504 individuals in the 1000 Genomes project phase 3. Now, we are only left to plot the result. There are many programs that can do this (like Excel) but I will show some R code (requiring ggplot2). If you are familiar with R, you can use the following script:
<pre>
# generate principal component plot
suppressPackageStartupMessages(library(ggplot2))
pca <- read.table('kgp.pca', header = TRUE)
pop <- read.table('kgp.pop', header = TRUE)
df <- merge(pca, pop)
pdf('kgp.pdf')
p <- list()
p[[1]] <- ggplot(df, aes(x=PC1, y=PC2))
p[[2]] <- ggplot(df, aes(x=PC1, y=PC3))
p[[3]] <- ggplot(df, aes(x=PC2, y=PC3))
p[[4]] <- ggplot(df, aes(x=PC1, y=PC4))
p[[5]] <- ggplot(df, aes(x=PC2, y=PC4))
p[[6]] <- ggplot(df, aes(x=PC3, y=PC4))
for (i in 1:6) {
print(p[[i]] + geom_point(aes(color=POP, shape=POP), alpha=1/3) +
scale_color_discrete() +
scale_shape_manual(values=0:(length(levels(df[,'POP']))-1)%%25+1) +
theme_bw())
}
dev.off()
</pre>
The R code will generate six images, the first being the following (see <a href="http://www.internationalgenome.org/category/population/">here</a> for a full legend):
<div class="separator" style="clear: both; text-align: center;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEioSu0GuXI4vAwrgY0xoOFrndOZisiKDYuRV2mdctM4an7XTl6BUBdxigt8bCOVf8JD1bldpjUsnWxARLB0om4U859uVjPhOTEAlGOC6mMaMZvsHU204fKsGOo08e65uxpkbO2bsAGSxZI/s1600/kgp.png" imageanchor="1" style="margin-left: 1em; margin-right: 1em;"><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEioSu0GuXI4vAwrgY0xoOFrndOZisiKDYuRV2mdctM4an7XTl6BUBdxigt8bCOVf8JD1bldpjUsnWxARLB0om4U859uVjPhOTEAlGOC6mMaMZvsHU204fKsGOo08e65uxpkbO2bsAGSxZI/s1600/kgp.png" /></a></div>
You might notice that the European populations on the top left, the African populations on the right, the East Asian populations in the bottom left, and the American and South Asian populations in the left, in between Europeans and East Asians. However, this plot does not make justice to American and South Asian populations.
A three-dimensional visualization would provide a much better understanding of how these five super-populations relate to each other. I will not share the ugly Matlab code I have used, but I will share the result:
<iframe width="560" height="315" src="https://youtube.com/embed/XkCCE8jC1F4?autoplay=1&loop=1&playlist=XkCCE8jC1F4" frameborder="0" allowfullscreen></iframe>
As you might now notice, the South Asian populations have their own non-overlapping location in space. You might also appreciate many other previously hiding details. The BEB population (Bengali from Bangladesh) now clearly leans towards the East Asian population, as historical accounts would expect. You will now notice very clearly other smaller details as the five ASW (African American) individuals with significant amount of Native American ancestry and one PEL (Peruvians from Lima) individual with significant amounts of East Asian ancestry. This one individual is <a href="http://catalog.coriell.org/0/Sections/Search/Sample_Detail.aspx?Ref=HG01944">HG01944</a> and he clearly deserves some extra attention. We will talk about him again in the next post. ;-)Unknownnoreply@blogger.com3tag:blogger.com,1999:blog-2735750853096182491.post-35256006272137305472015-10-24T20:50:00.003-07:002015-10-25T19:00:32.553-07:00Distinguishing parents from children from genotyping array dataFrom genotype data for a trio (father, mother, child), it is straightforward to identify whom the parents are. It suffices to identify which pair shares genotypes from the third person so that the Mendelian error rate is minimized. From genotype data for a duo (parent, child), it is straightforward to identify that the two individuals are immediately related and that the relationship is a parent-child one. However, determining who is the parent and who is the child is a challenge of a different nature.
<br />
<br />
Without phasing information, sharing between a parent and a child is quite a symmetric relationship. One piece of information that can give directionality would be the identification of de-novo mutations, or mutations not shared by an outlier individual, that would only be present in the child and not in the parent. From genotype data we cannot assay de-novo mutations, however.
<br />
<br />
The other piece of information that can give directionality is more closely related to linkage analysis. A distant relative might share some DNA with the parent and possibly the same amount of DNA with the child, or occasionally only part of that amount with the child. The reverse is also possible but quite unlikely.
<br />
<br />
If you have DNA information for a duo either in 23andMe or in AncestryDNA, you can test whether this information is actionable by first downloading your genotype data using the instructions from my previous <a href="http://apol1.blogspot.com/2015/09/visualize-dna-matches-graph.html">post</a> and then visualizing estimated sharing between members of the duo and shared distant relatives with an additional python <a href="https://github.com/freeseek/getmydnamatches">script</a>.
<br />
<br />
This is what I have obtained using my own 23andMe and AncestryDNA accounts:
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh51JM4KZrlVCLiLl1zbAmGxDLawHnN1vJ65AP02QyvAvxAb60FOH3FjR3y5THDbFZ96KDj6NGws-65rmACYsRTnOwiahjSPm3S6q4SInxAc42F7Yz2AQ9dr2dKu-_MLiGwcOvm9aoPQeI/s1600/share.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEh51JM4KZrlVCLiLl1zbAmGxDLawHnN1vJ65AP02QyvAvxAb60FOH3FjR3y5THDbFZ96KDj6NGws-65rmACYsRTnOwiahjSPm3S6q4SInxAc42F7Yz2AQ9dr2dKu-_MLiGwcOvm9aoPQeI/s1600/share.png" width=560 /></a>
Red markers indicate individuals for whom sharing differs with the individuals of the duo by an amount larger than 0.1% for 23andMe, or 5 centiMorgans for AncestryDNA. As you can see from the figure, there are always more red markers on the parent's side of the diagonal.
<br />
<br />
There are two explanations for this observation:<br />
1) A distant relative shares more than one segment with the parent and a smaller number of segments with the child<br />
2) A distant relative shares one large segment of DNA with the parent and a smaller chunk of that segment with the child due to a recombination event which split the segment in two parts, of which only one was passed to the child<br />
<br />
Due to the large number of distant cousins sharing only one segment of DNA, I believe that it is more common to observe case (2) although you still need to expect the passed segment to be at least 5 centiMorgans long in order to be detected and reported.
<br />
<br />
Notice also that occasionally the child might share more with the distant cousin than the parent does (red markers under the diagonal). It is possible that statistical noise causes a shared DNA segment to be reported in the child but not in the parent. It is also occasionally possible that the distant relative is related to the child both through the mother and through the father, though in my experience this circumstance seems to be rare.
<br />
<br />
<b>Limitations of 23andMe</b>
<br />
<br />
While the Relative Finder section of 23andMe shows a long list of DNA matches for each profile, the script used in this post will not pair matches showing up as anonymous. This most likely will exclude 60-70% DNA matches, depending on profile anonymity chosen by your distant cousins.Unknownnoreply@blogger.com4tag:blogger.com,1999:blog-2735750853096182491.post-10619427141413319252015-09-27T11:43:00.000-07:002015-10-05T14:43:50.331-07:00Visualize DNA matches graphThrough high-throughput SNP genotyping it has become possible to identify shared DNA with our fifth degree cousins. Two companies, <a href="https://23andme.com/">23andMe</a> and <a href="http://dna.ancestry.com/">AncestryDNA</a>, currently each hold each a database of more than one million people of dense genotype data. Each database gives us the possibility to identify hundreds of distant DNA cousins as each service reports a list of all individuals who share DNA segments at least 5 centimorgans long.
<br />
<br />
DNA cousins from 1st up to 4th degree will usually share several such DNA segments making them stand out from the crowd of more distantly related DNA cousins. For the large list of individuals that only share with us one or two DNA segments, it is difficult to distinguish between a match as close as a fifth degree cousin and someone much more distantly related who just happens to share a long DNA segment which through luck stayed intact across several centuries.
<br />
<br />
Is there any information other than the number of shared DNA segments and their length in centimorgans that we can leverage to prioritize our DNA matches? Since <a href="http://blogs.ancestry.com/ancestry/2015/08/26/see-your-dna-matches-in-a-whole-new-way/">August 26th 2015</a>, AncestryDNA has launched a feature whereby for a given match we can further learn whether it in turn shares DNA with other of our DNA matches, with some limitations. We can also do the same for DNA matches that have agreed to share genomes with us on 23andMe, though it is a little bit more cumbersome.
<br />
<br />
We can therefore define a <i>DNA matches graph</i> where each node is one of our DNA matches and two nodes are connected by an edge if the two corresponding DNA matches share DNA with each other. What can we learn from this graph? What does the topology look like? Does the graph contain connected components? For example, if our parents were to be from very different places, it would be quite unlikely that our DNA matches from our mother's side would be connected to our DNA matches from our father's side. The same holds if our grandparents are from different places, and so on.
<br />
<br />
To be able to glance at our <i>DNA matches graph</i>, we need to access all the information that either AncestryDNA or 23andMe provides us at once. To retrieve this information manually is just impractical. To overcome this hurdle, I have put together a few python <a href="https://github.com/freeseek/getmydnamatches">scripts</a> which will do the job for us.
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2O6uzx_GibIsaTRYrkTJP6ho_n68zMh6eOQ7kG40t0eE9vR6lXZM8Q9_Pz8hGXfj9zlMaNS-7gI2hyphenhyphen_rLVXAEtftlC8BZXEbR24Sqbff4Tec9qq7h6_mBL8k1X0Y3nYLOLuL3Z7GZ830/s1600/getmydna.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj2O6uzx_GibIsaTRYrkTJP6ho_n68zMh6eOQ7kG40t0eE9vR6lXZM8Q9_Pz8hGXfj9zlMaNS-7gI2hyphenhyphen_rLVXAEtftlC8BZXEbR24Sqbff4Tec9qq7h6_mBL8k1X0Y3nYLOLuL3Z7GZ830/s1600/getmydna.png" /></a>
<br />
The workflow runs as follows:<br />
1) We dump all DNA matches information in our user account
<br />
2) We process the information to generate a table of all pairwise DNA sharing
<br />
3) We visualize the <i>DNA matches graph</i>
<br />
<br />
However, for now, you will have to be familiar with running python commands and installing a few necessary python packages. Here are some examples of what came out for my mother's 23andMe test:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEigsm-44HMIqDMkN5vzYMlM-5QR8R3jVrC2Au4GIbxm3Sk_UVH7PjxJfPLVl6wfGQ1-tR_jqMI7QLnfLQqWnuDBAkcFy2g6awq4-scuDXwys2c7vqzN58c0qhEF4jSzLKbjzOx2yr0sA6E/s1600/mother.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEigsm-44HMIqDMkN5vzYMlM-5QR8R3jVrC2Au4GIbxm3Sk_UVH7PjxJfPLVl6wfGQ1-tR_jqMI7QLnfLQqWnuDBAkcFy2g6awq4-scuDXwys2c7vqzN58c0qhEF4jSzLKbjzOx2yr0sA6E/s1600/mother.png" /></a>
<br />
The blue nodes correspond to DNA matches on my paternal grandmother's side and the pink nodes correspond to DNA matches on my maternal grandmother's side. Interestingly these two groups form two very separate clusters and this is to be expected as my maternal grandparents come from two Italian towns quite distant from each other.
<br />
<br />
The following example is from the AncestryDNA account of an American friend of mine:
<br />
<a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSJsKiWlNxpN51eNmYmJKzVkvw2ZjVqaJkOK4O47mCsbnBadwj45iPg0_0bQBr_GO4BaCYMr2G4Bjvng4Lc64NR3cFekhoqfPUX_6NGdCK6tZk9bPi5JfF-1mFS06WRmtVNJGMmKPY4Hs/s1600/friend.png" imageanchor="1" ><img border="0" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEjSJsKiWlNxpN51eNmYmJKzVkvw2ZjVqaJkOK4O47mCsbnBadwj45iPg0_0bQBr_GO4BaCYMr2G4Bjvng4Lc64NR3cFekhoqfPUX_6NGdCK6tZk9bPi5JfF-1mFS06WRmtVNJGMmKPY4Hs/s1600/friend.png" /></a>
<br />
The small nodes correspond to DNA matches estimated as my friend's distant cousins, while the larger nodes are estimated as 4th cousins or closer, and the dark colored nodes correspond to nodes for which AncestryDNA has identified the most recent common ancestor in the family trees within nine generations. The clusters make some sense. For example, we were able to identify a little one corresponding to a little town from Poland where my friend's great-great grandmother is from. It is a great graph and I am jealous because if you have Northern European roots, you get a lot more DNA matches than if you have Southern European roots like me, most likely due to the bias in which individuals buy the AncestryDNA test.
<br />
<br />
If you are savvy enough, I encourage you to use the tool on your own data and do send me feedback! :-)
<br />
<br />
<b>Limitations of 23andMe and AncestryDNA</b>
<br />
<br />
When taking a DNA test with 23andMe or AncestryDNA all matches sharing at least 5 centimorgans with the individual being tested will be reported for further review. In either case some matches will be anonymous and some matches will share as much as their own family tree, depending on the settings selected by the user. I believe for chromosome X matches for males 23andMe might be lowering the threshold to something like 2 centimorgans as it is significantly easier to identify sharing on the X chromosome among males. This 5 centimorgans threshold for shared segments is likely selected for two reasons:
<br />
<br />
1) Shorter segments become increasingly likely to be false positives below 5 centimorgans without your genome being properly <i>phased</i> which is mostly not the case with the exception of those customers who also had their parents tested.
<br />
<br />
2) Shorter segments are likely to be due to really far away relationships and as such they might have very little value for genealogy as they would trace back several centuries, much beyond what the paper trail allows in most countries.
<br />
<br />
<b>Limitations of 23andMe</b>
<br />
<br />
When browsing the <i>Relative Finder</i> feature within 23andMe, you can see a full list of individuals sharing with you DNA segments larger than 5 centimorgans. For each DNA match you can also see how many separate segments you share, unlike AncestryDNA. However, you cannot determine whether two of your DNA matches are DNA matches with each other unless both users have agreed to share their genetic profiles with you.
<br />
<br />
Due to this limitation you can only build the <i>DNA matches graph</i> for a restricted class of individuals who are engaged enough in the service to be able to manage their sharing requests. Furthermore, you can only retrieve information for one pair of individuals at a time, making the process completely impractical if attempted manually. Due to how slow the 23andMe server is at retrieving such information, even with the script mentioned above, it will take about three seconds to retrieve information about each pair. For 250 individuals it might take approximately 24 hours, even if most retrieval requests will yield empty results.
<br />
<br />
<b>Limitations of AncestryDNA</b>
<br />
<br />
Unlike 23andMe, AncestryDNA does not report the number of shared segments among DNA matches nor the amount of genome shared. However, when browsing the site the browser does receive two variables invisible to the user:
<br />
<br />
1) <i>sharedCentimorgans</i> being an estimate of the amount of DNA shared across shared DNA segments larger than 5 centimorgans
<br />
<br />
2) <i>meiosisValue</i> being an estimate from <i>sharedCentimorgans</i> of the number of meioses separating you and your DNA match and this number is never larger than 10
<br />
<br />
From what I observed pairs of individuals within the AncestryDNA database can then be split in three categories:
<br />
<br />
1) Pairs that share no DNA segments larger than 5 centimorgans
<br />
<br />
2) Pairs that share at least one DNA segment longer than 5 centimorgans with a <i>meiosisValue</i> of exactly 10 (corresponding to <i>sharedCentimorgans</i> between 5 and 17.65, that is, relationships classified as 10:<i>DISTANT COUSIN</i>)
<br />
<br />
3) Pairs that share at least one DNA segment longer than 5 centimorgans with a <i>meiosisValue</i> smaller than 10 (corresponding to <i>sharedCentimorgans</i> greater than 17.65, that is, relationships classified as 0:<i>SELF/TWIN</i>, 1:<i>PARENT/CHILD</i>, 2:<i>IMMEDIATE FAMILY</i>, 3:<i>CLOSE FAMILY</i>, 4:<i>1ST COUSIN</i>, 5,6:<i>2ND COUSIN</i>, 7:<i>3RD COUSIN</i>, 8,9:<i>4TH COUSIN</i>)
<br />
<br />
When accessing your AncestryDNA test, you will have reported all DNA matches corresponding to the last two categories. Most (in my case 96-98%) of your pairs including you and one of your DNA matches will fall in the second category. However, when you are reviewing one of your DNA matches <i>A</i> in AncestryDNA, the Shared Matches feature will list only DNA matches <i>B</i> such that <b>both</b> <i>A</i> and <i>B</i> <b>and</b> you and <i>B</i> fall in the third category as also explained in another <a href="http://www.thegeneticgenealogist.com/2015/08/28/ancestrydna-announces-new-in-common-with-tool/">post</a>.
<br />
<br />
As of this writing, even when listing the DNA matches you have in common with your father/mother you are actually only listing those DNA matches which make a third category pair with your father/mother. Such limitation is not present with 23andMe but I believe this is just a bug on AncestryDNA's side which will get fixed soon.
<br />
<br />
To summarize, through AncestryDNA you will only have access to those edges of the <i>DNA matches graph</i> corresponding to pairs sharing at least 17.65 centimorgans of DNA. This is quite a conservative selection and likely the result of careful debate within AncestryDNA to find a balance between being at the same time informative and concise.Unknownnoreply@blogger.com3tag:blogger.com,1999:blog-2735750853096182491.post-65049709610102911542015-05-19T12:35:00.000-07:002016-01-23T15:25:24.073-08:00Annotate CpG mutations in a VCF fileMutations in the human genome are not made all the same. Even when restricting our attention to the most common form of human variation, that is, single nucleotide polymorphism, there are different categories. On a first approximation we have transversion and transitions. And among transitions we have transitions on CpG sites, which are much more common than transitions on non-CpG sites. Hence, the need to distinguish the two. If you have mutations encoded in the standard VCF format, it is easy to distinguish transitions from transversions. The first kind shows as a A<->G or a C<->T, while the second kind shows as a A<->C, A<->T, C<->G, or G<->T. But distinguishing CpG transitions from non-CpG transitions within a VCF is not possible without additional information as we need information about the sequence context.<br />
<br />
We definitely need the reference genome sequence of the organism the VCF file refers to. However there is no tool to my knowledge that will take a VCF file and a fasta sequence and yield as output an annotated VCF file indicating which mutations are CpG mutations (it would be nice if someone wrote a bcftools pluging for this purporse). But there are tools that will annotate a VCF file given the presence of a mutation in another VCF file. Therefore one solution would be to have a VCF file with all possible CpG mutations. Here is how to do that:<br />
<pre><code>
# create your personal binary directory
mkdir -p ~/bin/
# download human genome reference (GRCh37)
wget -O- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz | gunzip > human_g1k_v37.fasta &&
samtools faidx human_g1k_v37.fasta
# install latest version of bedtools
git clone https://github.com/arq5x/bedtools2.git &&
cd bedtools2 && make && cd .. &&
cp bedtools2/bin/bedtools ~/bin/
# create a VCF file with all CpG mutations (it takes several hours)
awk '{for (i=1; i<$2; i++)
print $1"\t"i-1"\t"i+1"\t"$1" "i}' human_g1k_v37.fasta.fai |
bedtools getfasta -name -fi human_g1k_v37.fasta \
-bed /dev/stdin -fo /dev/stdout |
tr '\n' ' ' | tr '>' '\n' | grep "CG $" |
awk -v OFS="\t" 'BEGIN {print "##fileformat=VCFv4.1";
print "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO"}
{print $1"\t"$2"\t.\tC\tT\t.\t.\t.";
print $1"\t"$2+1"\t.\tG\tA\t.\t.\t."}' |
bgzip > cpg.vcf.gz &&
tabix -f cpg.vcf.gz
</code></pre>
To be completely thorough these mutations will be CpG mutations assuming that the reference sequence is the ancestral sequence which is not always the case. But for the majority of rare mutations it is a fair assumption. Furthermore, this code does not take into account whether a CpG mutation is located within a CpG island and this is also information that might be important as not all CpG sites are equally mutable.<br />
<br />
Now that we have our VCF file with all CpG mutations, we can use SnpSift to annotate our input VCF file. I selected SnpSift because it is quite fast and flexible compared to other tools like bcftools (see <a href="https://github.com/samtools/bcftools/issues/108">here</a>). This can be achieved as follows:<br />
<pre><code>
# install latest version of snpEff/SnpSift
wget http://downloads.sourceforge.net/project/snpeff/snpEff_latest_core.zip &&
unzip snpEff_latest_core.zip
# annotate your VCF file with SnpSift
java -jar snpEff/SnpSift.jar annotate -exists CPG cpg.vcf.gz input.vcf.gz | bgzip > output.vcf.gz &&
tabix -f output.vcf.gz
</code></pre>
If you want to extract from the VCF file only the variants which are CpG mutations, the following code will work:<br />
<pre><code>
# install latest version of bcftools
git clone --branch=develop git://github.com/samtools/bcftools.git
git clone --branch=develop git://github.com/samtools/htslib.git
cd htslib && make && cd ..
cd bcftools && make && cd ..
mv bcftools/bcftools ~/bin/
# extract from VCF file all CpG mutations
bcftools view -Oz -i "CPG==1" output.vcf.gz -o output.cpg.vcf.gz &&
tabix -f output.cpg.vcf.gz
</code></pre>Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-2735750853096182491.post-3682322057956173892014-11-25T12:02:00.000-08:002015-07-07T13:40:32.115-07:00Best practice for converting VCF files to plink formatConverting VCF files to plink format has never been easier. However, there are a few issues related to some intrinsic limitations of the plink format. The first is related to the fact that variants in a plink file are bi-allelic only, while variants in a VCF file can be multi-allelic. The second is related to an intrinsic limitation of plink which makes indel definitions ambiguous. Here is an example: is the following variant an insertion or a deletion compared to the GRCh37 reference?<br />
<pre><code>
20 31022441</code> A AG</pre>
<br />
There is no way to tell, as the plink format does not record this information.<br />
<br />
Keeping this in mind, we are going to need two pieces of software for the conversion, bcftools and plink2. Here how to install them:<br />
<pre><code>
# create your personal binary directory
mkdir -p ~/bin/
# install latest version of bcftools
git clone --branch=develop git://github.com/samtools/bcftools.git
git clone --branch=develop git://github.com/samtools/htslib.git
cd htslib && make && cd ..
cd bcftools && make && cd ..
mv bcftools/bcftools ~/bin/
# install latest version of plink2
wget https://www.cog-genomics.org/static/bin/plink/plink_linux_x86_64.zip
unzip -od ~/bin/ plink_linux_x86_64.zip plink
</code></pre>
<br />
We are also going to need a copy of the GRCh37 reference. If you don't have this, it will take a while to download, but it can be done with the following commands:<br />
<pre><code>
wget -O- ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/human_g1k_v37.fasta.gz | gunzip > human_g1k_v37.fasta
samtools faidx human_g1k_v37.fasta
</code></pre>
<br />
The following command will take the VCF file, strip the variant IDs, split multi-allelic sites into bi-allelic sites, assign names to make sure indels will not become ambiguous, and finally convert to plink format:<br />
<pre><code>
bcftools norm -Ou -m -any input.vcf.gz |
bcftools norm -Ou -f human_g1k_v37.fasta |
bcftools annotate -Ob -x ID \
-I +'%CHROM:%POS:%REF:%ALT' |
plink --bcf /dev/stdin \
--keep-allele-order \
--vcf-idspace-to _ \
--const-fid \
--allow-extra-chr 0 \
--split-x b37 no-fail \
--make-bed \
--out output
</code></pre>
<br />
Usually dbSNP IDs are used for variants in VCF files, but in my opinion this is really a bad practice especially after splitting multi-allelic sites into bi-allelic sites (see <a href="http://annovar.openbioinformatics.org/en/latest/articles/dbSNP/">here</a> for a discussion).<br />
<br />
The first command will split multi-allelic alleles so that a variant like the following:<br />
<pre><code>
20 31022441 AG AGG,A
</code></pre>
<br />
Will become two variants as follows:<br />
<pre><code>
20 31022441 AG AGG
20 31022441 AG A
</code></pre>
<br />
The second command will make sure that after having been split, indels are normalized so to have a unique representation. The above multi-allelic variant would then become:<br />
<pre><code>
20 31022441 A AG
20 31022441 AG A
</code></pre>
<br />
Notice that plink will have a very hard time to distinguish between the two variants above, as they look alike once you forget which allele is the reference allele. The third command will assign a unique name to each bi-allelic variant:<br />
<pre><code>
20 31022441 20:31022441:A:AG A AG
20 31022441 20:31022441:AG:A AG A
</code></pre>
<br />
The fourth command will convert the final output to plink binary format. At this point, you will have a plink binary file with unique names for each variant. You can test that this is the case with the following command:<br />
<pre><code>
cut -f2 output.bim | sort | uniq -c | awk '$1>1'
</code></pre>
<br />
If the above command returns any output, then you still have duplicate IDs but this means your VCF file must be flawed to begin with. Plink will not like having multiple variants with the same ID.<br />
<br />
You should now be able to use and to merge the file you have generated with other plink files generated from VCF files in the same way.<br />
<br />
Your last worry is individuals' sex, as the VCF format, contrary to plink format, does not encode this information. If your VCF file contains enough information about the X chromosome, you should be able to assign the sex straight from genotype. This is a command that might do the trick:<br />
<br />
<code></code><br />
<pre><code>plink --bfile output --impute-sex .3 .3 --make-bed --out output && rm output.{bed,bim,fam}~
</code></pre>
<br />
However, the two thresholds you should use (.3 .3 in the above example) to separate males from females are going to be dependent on the particular type of data you are using (e.g. exome or whole genome) and might have to be selected ad hoc. See <a href="https://www.cog-genomics.org/plink2/basic_stats#check_sex">here</a> for additional information.
<br />
A python wrapper to perform the above operations is available online (<a href="https://raw.githubusercontent.com/freeseek/gpipe/master/vcf2plink.py">here</a>).Unknownnoreply@blogger.com21tag:blogger.com,1999:blog-2735750853096182491.post-35903555622707927452014-05-13T12:57:00.002-07:002014-05-13T12:57:32.522-07:00How do I identify all the homopolymer run indels in a VCF file?With sequencing projects we analysts usually are left to deal with genotype data in VCF format. These files can contain genotypes for a variety of markers, from SNPs to structural variants. Small structural variants, commonly known as indels, are usually harder to call from sequencing data, and we tend to be wary of their genotypes. However, not all indels are born the same. Among indels there are a variety of different kind of indels, and depending on your project, you might be more or less interested in these. A form of particularly difficult indels to call are homopolymer run indels, that is, deletion or insertions within repeats of the same nucleotide. A very famous homopolymer run indel is the one causing medullary cystic kidney disease type 1 by disrupting the coding sequence within a VNTR of MUC1.<br />
<br />
It might be important to identify in your VCF file which indel are homopolymer run indels, but this is not so straightforward. Also, this information cannot be extracted from the VCF file alone without the reference, so the information from two files needs to be integrated. I thought of a little hack to do this hastily without writing any serious amount of code using bedtools only. The following bash script should do it:<br />
<blockquote class="tr_bq">
#!/bin/bash<br />input="$1" # input VCF file in any format accepted by bcftools<br />ref="$2" # reference genome in fasta format </blockquote>
<blockquote class="tr_bq">
bcftools view -HG $input | cut -f1-5 | awk '$4!~"^[ACGT]($|,)" || $5!~"^[ACGT]($|,)"' |<br /> awk '{print $1"\t"$2+length($4)-1"\t"$2+length($4)+5"\t"$3}' |<br /> bedtools getfasta -name -fi $ref -bed /dev/stdin -fo /dev/stdout |<br /> awk 'NR%2==1 {printf substr($0,2)"\t"} NR%2==0' |<br /> grep "AAAAAA$\|CCCCCC$\|GGGGGG$\|TTTTTT$" | cut -f1 | sort</blockquote>
The script will take the VCF file and the reference as output and will yield the names of all markers identified as homopolymer run indels. Notice that it requires each marker in your VCF file to have its own unique ID. Also, the definition of homopolymer run indel is arbitrary here, consisting of the repeat of six base pairs. It should be easy to change the code if you want to tweak that.Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-2735750853096182491.post-40490893332559222772013-12-16T10:24:00.001-08:002013-12-16T10:24:23.517-08:00Estimating global ancestral components with 1000 GenomesFollowing the last two posts, here I will show you how to use the PCA results after merging your dataset with the 1000 Genomes project dataset to estimate global ancestral components for the four main ancestries around the world, that is, European, African, Native American, and Asian.<br />
<br />
Since there are no Native Americans in the 1000 Genomes project, we need the global ancestral estimates for the Latinos in the project which will be necessary to estimate this component. These estimates are available as part of the project, so we just need to download them:<br />
<br />
<pre>anc="%DIRECTORY WITH ANCESTRAL COMPONENTS%"
mkdir -p $anc
# download global proportions table
(echo "IID EUR AFR NAT UNK";
for pop in ASW CLM MXL PUR; do
wget -O- ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/phase1/analysis_results/ancestry_deconvolution/${pop}_phase1_globalproportions_and_unknown.txt
done) > $anc/gap
</pre>
<br />
The 1000 Genomes project individuals can then be used to train a linear model that estimates global ancestral components. If you followed the directory structure for the PCA as discussed in the previous <a href="http://apol1.blogspot.com/2013/12/pca-of-your-plink-dataset-with-1000.html">post</a>, then the following R script will take the three tables $kgp/$set.kgp.pop kgp/$set.kgp.pca $anc/gap as input and will output a fourth table, say kgp/$set.anc, with the global ancestry estimates.
<br />
<br />
<pre>#!/usr/bin/Rscript
args <- commandArgs(trailingOnly = TRUE)
</pre>
<pre># load and merge input tables
df1 <- read.table(args[1],header=TRUE)
df2 <- read.table(args[2],header=TRUE)
df3 <- read.table(args[3],header=TRUE)
data <- merge(df1,df2,all=TRUE,incomparables=NaN)
data <- merge(data,df3,all=TRUE,incomparables=NaN)
data$ASN <- NaN
# European populations
eur <- data$POP=="CEU" | data$POP=="FIN" | data$POP=="GBR" | data$POP=="IBS" | data$POP=="TSI"
data[which(eur),"EUR"] <- 1; data[which(eur),"AFR"] <- 0; data[which(eur),"NAT"] <- 0; data[which(eur),"ASN"] <- 0
# African populations
afr <- data$POP=="YRI"
data[which(afr),"EUR"] <- 0; data[which(afr),"AFR"] <- 1; data[which(afr),"NAT"] <- 0; data[which(afr),"ASN"] <- 0
# Asian populations
asn <- data$POP=="CDX" | data$POP=="CHB" | data$POP=="CHD" | data$POP=="CHS" | data$POP=="JPT"
data[which(asn),"EUR"] <- 0; data[which(asn),"AFR"] <- 0; data[which(asn),"NAT"] <- 0; data[which(asn),"ASN"] <- 1
# Admixed populations
mix <- data$POP=="ASW" | data$POP=="CLM" | data$POP=="MXL" | data$POP=="PUR"
for (anc in c("EUR","AFR","ASN")) {
data[which(mix),anc] <- data[which(mix),anc] * 1/(1-data[which(mix),"UNK"])
}
data[which(mix),"ASN"] <- 0
# African American samples to be excluded
exclude = c("NA20299","NA20314","NA20414","NA19625","NA19921")
for (anc in c("EUR","AFR","NAT","ASN")) {
data[is.element(data$IID,exclude),anc] <- NaN
}
# predict global ancestral proportions
for (anc in c("EUR","AFR","NAT","ASN")) {
tmp <- lm(as.formula(paste0(anc," ~ PCA1 + PCA2 + PCA3")), data)
data[is.na(data[,anc]),anc] = predict.lm(tmp,data[is.na(data[,anc]),])
}
# output data to disk
write.table(data[,c("FID","IID","EUR","AFR","NAT","ASN")],args[4],row.names=FALSE,quote=FALSE)
</pre>
<br />
If you look into the code, you will see that some African American samples are excluded from the computation. This is due to the fact that a few African American individuals in the 1000 Genomes project have a significant amount of Native American ancestry, but this does not show up in the global ancestry proportion tables as it was not modeled. It is therefore better to exclude those individuals from these computations.<br />
<br />
Also, you might notice that, due to the inherent noise in the way the global ancestral components are modeled, it is possible that some ancestral components will be negative. This is to be expected, as these are noisy estimates. Make sure any further downstream analysis takes care of this eventuality. Overall, estimating ancestral components from principal components is a more noisy business than estimating ancestral components from local ancestry deconvolutions. However, for many purposes these estimates should work just fine.Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-2735750853096182491.post-76645647388681513282013-12-13T13:14:00.001-08:002013-12-13T13:14:19.118-08:00PCA of your plink dataset with 1000 GenomesEstimating global ancestral components of your samples should really not be harder than running a principal component analysis (PCA). Here in this post, I will show you how you can do exactly that, using the 1000 Genomes project samples. You will first need your dataset and the 1000 Genomes dataset in plink format, as explained in the <a href="http://apol1.blogspot.com/2013/12/1000-genomes-pca-analysis.html">previous post</a><br />
<br />
As usual, you will need a few variables defined:
<br />
<pre>set="%PLINK DATASET PREFIX%"
pth="%PATH TO THE 1000 GENOMES PLINK DATASET%"
gct="%DIRECTORY WITH GCTA%"
</pre>
<br />
It is a good idea to make sure that you dataset is properly cleaned and pruned before merging the data with the 1000 Genomes project data. It is also important that SNPs be coded in the forward strand, as this is how the 1000 Genomes project is coded as. The following code will identify SNPs that are shared and compatible between the datasets, merge the individuals over these SNPs, and compute the genetic relationship matrix necessary for the PCA analysis:
<br />
<br />
<pre>mkdir -p kgp
# identify SNPs compatible with the merge
awk 'NR==FNR {chr[$2]=$1; pos[$2]=$4; a1[$2]=$5; a2[$2]=$6} NR>FNR && ($2 in chr) && chr[$2]==$1 && pos[$2]==$4 && (a1[$2]==$5 && a2[$2]==$6 || a1[$2]==$6 && a2[$2]==$5) {print $2}' $set.bim <(cat $(seq -f $pth/kgp.chr%g.bim 1 22) $pth/kgp.chrX.bim) > kgp/extract.txt
# extract SNPs from the 1000 Genomes dataset
for chr in {1..22} X; do plink --noweb --bfile $pth/kgp.chr$chr --extract kgp/extract.txt --make-bed --out kgp/kgp.chr$chr; done
/bin/cp kgp/kgp.chr1.fam kgp/kgp.fam
for chr in {1..22} X; do cat kgp/kgp.chr$chr.bim; done > kgp/kgp.bim
(echo -en "\x6C\x1B\x01"; for chr in {1..22} X; do tail -c +4 kgp/kgp.chr$chr.bed; done) > kgp/kgp.bed
for chr in {1..22} X; do /bin/rm kgp/kgp.chr$chr.bed kgp/kgp.chr$chr.bim kgp/kgp.chr$chr.fam kgp/kgp.chr$chr.log; done
# merge two datasets
plink --noweb --bfile $set --bmerge kgp/kgp.bed kgp/kgp.bim kgp/kgp.fam --make-bed --out kgp/$set.kgp
# create population categorical data
awk 'BEGIN {print "FID","IID","POP"} NR==FNR {pop[$2]=$7} !($2 in pop) {pop[$2]="SET"} NR>FNR {print $1,$2,pop[$2]}' $pth/20130606_g1k.ped kgp/$set.kgp.fam > kgp/$set.kgp.pop
# use GCTA to compute the genetic relationship matrix
gcta64 --bfile kgp/$set.kgp --autosome --make-grm-bin --out kgp/$set.kgp --thread-num 10
</pre>
At this point we are ready to compute principal components.
<br />
<pre># compute principal components using GCTA
gcta64 --grm-bin kgp/$set.kgp --pca 20 --out kgp/$set.kgp --thread-num 10
# create PCA matrix for use with R data.frame
(echo "FID IID "`seq -f PCA%.0f -s" " 1 20`;
(echo -n "#eigvals: "; head -n20 kgp/$set.kgp.eigenval | tr '\n' ' '; echo; cat kgp/$set.kgp.eigenvec) |
awk 'NR==1 {n=NF-1; for (i=1; i<=n; i++) eval[i]=$(i+1)}
NR>1 {NF=2+n; for (i=1; i<=n; i++) $(i+2)*=sqrt(eval[i]); print}') > kgp/$set.kgp.pca
</pre>
<br />
The previous code should have created two files, kgp/$set.kgp.pop and kgp/$set.kgp.pca, which can be easily loaded into R with the following script which takes the two files as paramters:
<br />
<br />
<pre>#!/usr/bin/Rscript
suppressPackageStartupMessages(library(ggplot2))
args <- commandArgs(trailingOnly = TRUE)
df1 <- read.table(args[1], header = TRUE)
df2 <- read.table(args[2], header = TRUE)
data <- merge(df1,df2)
ggplot(data,aes(x=PCA1,y=PCA2,colour=POP,shape=POP)) + geom_point(alpha=1/3) + scale_colour_discrete() + scale_shape_manual(values=1:length(levels(data[,"POP"])))
</pre>
<br />
Notice that you will need R and ggplot2 already installed to run this script.Unknownnoreply@blogger.com4tag:blogger.com,1999:blog-2735750853096182491.post-89868005232342206312013-12-12T16:17:00.000-08:002014-10-29T13:08:25.796-07:001000 Genomes PCA analysisThe easiest way run a PCA analysis with the 1000 Genomes samples is to download the data, convert it to plink format, and use GCTA to perform the bulk of the computation. The dataset provided on the beagle website is likely the easiest to start with. Here some code that can get the work done.<br />
<br />
First we need to setup a few variables:
<br />
<pre>gct="%DIRECTORY WITH GCTA%"
vcf="%DIRECTORY WITH KGP VCF FILES%"
</pre>
<br />
And then we can download all necessary files:
<br />
<pre># install plink and utilities to convert to plink format
sudo apt-get install plink
wget http://broadinstitute.org/~giulio/vcf2plink/vcf2tfam.sh
wget http://broadinstitute.org/~giulio/vcf2plink/vcf2tped.sh
chmod a+x vcf2tfam.sh vcf2tped.sh</pre>
<pre># download the 1000 Genomes files
pfx="http://bochet.gcc.biostat.washington.edu/beagle/1000_Genomes_phase1_vcf/chr"
sfx=".1kg.ref.phase1_release_v3.20101123.vcf.gz"
for chr in {1..22} X; do wget $pfx$chr$sfx $pfx$chr$sfx.tbi -P $vcf/; done</pre>
<pre>wget ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/working/20130606_sample_info/20130606_g1k.ped
# install GCTA
wget http://www.complextraitgenomics.com/software/gcta/gcta_1.22.zip
mkdir -p $gct
unzip gcta_1.22.zip gcta64 -d $gct/
</pre>
<br />
For each chromosome we now need to create a plink binary file:
<br />
<pre># convert each VCF file to a plink binary file with sex information included
pfx="$vcf/chr"
for chr in {1..22} X; do
tabix -H $pfx$chr$sfx | ./vcf2tfam.sh /dev/stdin > kgp.chr$chr.tfam
tabix $pfx$chr$sfx $chr | ./vcf2tped.sh /dev/stdin > kgp.chr$chr.tped
p-link --noweb --tfile kgp.chr$chr --make-bed --out kgp.chr$chr
/bin/rm kgp.chr$chr.tfam kgp.chr$chr.tped kgp.chr$chr.nosex kgp.chr$chr.log
awk 'NR==FNR {sex[$2]=$5} NR>FNR {$5=sex[$2]; print}' 20130606_g1k.ped kgp.chr$chr.fam > kgp.chr$chr.fam.tmp && /bin/mv kgp.chr$chr.fam.tmp kgp.chr$chr.fam
done
</pre>
<br />
At this point we can use plink to prune each chromosome and then join the individual datasets into a single large dataset as discussed in a <a href="http://apol1.blogspot.com/2013/10/merging-plink-binary-files.html">previous post</a>:<br />
<pre># prune each chromosome
for chr in {1..22} X; do
p-link --noweb --bfile kgp.chr$chr --maf .01 --indep 50 5 2 --out kgp.chr$chr
p-link --noweb --bfile kgp.chr$chr --extract kgp.chr$chr.prune.in --make-bed --out kgp.chr$chr.prune
done
# join pruned datasets into single dataset
cat kgp.prune.chr1.fam > kgp.prune.fam
for chr in {1..22} X; do cat kgp.chr$chr.prune.bim; done > kgp.prune.bim
(echo -en "\x6C\x1B\x01"; for chr in {1..22} X; do tail -c +4 kgp.chr$chr.prune.bed; done) > kgp.prune.bed
</pre>
<br />
We can now run GCTA to compute the genetic relationship matrix and then the main principal components:<br />
<pre># perform PCA analysis of the whole dataset
$gct/gcta64 --bfile kgp.prune --autosome --make-grm-bin --out kgp --thread-num 10
$gct/gcta64 --grm-bin kgp --pca 20 --out kgp --thread-num 10
# create PCA matrix to be loaded into R
(echo "FID IID "`seq -f PCA%.0f -s" " 1 20`;
(echo -n "#eigvals: "; head -n20 kgp.eigenval | tr '\n' ' '; echo; cat kgp.eigenvec) |
awk 'NR==1 {n=NF-1; for (i=1; i<=n; i++) eval[i]=$(i+1)}
NR>1 {NF=2+n; for (i=1; i<=n; i++) $(i+2)*=sqrt(eval[i]); print}') > kgp.pca
</pre>
<br />
This will generate a table that can be easily loaded into R for further analyses.Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-2735750853096182491.post-7097899305596368532013-11-18T14:24:00.001-08:002013-12-02T11:11:40.077-08:00IBD pipelineThere are a lot of interesting analyses that segments inherited identical-by-descent (IBD) can be informative for. However, it takes some steps to setup a fully functional IBD pipeline. Here is a basic one that will work with any dataset in binary plink format. It will first phase the data with hapi-ur, then detect IBD segments using beagle, and finally detect clusters of IBD individuals using dash.<br />
<br />
To begin with, you need to setup a few variables:
<br />
<pre>hpr="%DIRECTORY WITH HAPI-UR BINARY FILES%"
gmp="%DIRECTORY WITH GENETIC MAP FILES%"
bgl="%DIRECTORY WITH BEAGLE BINARY FILE%"
dsh="%DIRECTORY WITH DASH BINARY FILE%"
set="%PLINK DATASET PREFIX%"
</pre>
<br />
And have a couple of programs pre-installed on your Linux machine:
<br />
<pre># install plink
sudo apt-get install plink
# install hapi-ur
wget https://hapi-ur.googlecode.com/files/hapi-ur-1.01.tgz
tar xzvf hapi-ur-1.01.tgz
mkdir -p $hpr
/bin/cp hapi-ur-1.01/hapi-ur hapi-ur-1.01/insert-map.pl $hpr/
# install genetic maps
wget http://mathgen.stats.ox.ac.uk/impute/genetic_maps_b37.tgz
tar xzvf genetic_maps_b37.tgz
mkdir -p $gmp
for i in {1..22} X_PAR1 X_nonPAR X_PAR2; do gzip -c genetic_maps_b37/genetic_map_chr${i}_combined_b37.txt > $gmp/chr$i.gmap.gz; done
# install beagle
mkdir -p $bgl
wget http://faculty.washington.edu/browning/beagle/b4.r1196.jar -P $bgl/
# install dash
wget http://www1.cs.columbia.edu/~gusev/dash/dash-1-1.tar.gz
tar xzvf dash-1-1.tar.gz
mkdir -p $dsh
/bin/cp dash-1-1/bin/dash_cc $dsh/
</pre>
<br />
The following pipeline will run the analysis independently on each autosome:
<br />
<pre># perform IBD and DASH analysis separately for each autosome
for i in {1..22}; do
# create a directory
mkdir -p chr$i
# extract autosome of interest
p-link --noweb --bfile $set --chr $i --make-bed --out chr$i/$set
# generate marker file with genetic distance
zcat $gmp/chr$i.gmap.gz | $hpr/insert-map.pl chr$i/$set.bim - | awk -v OFS="\t" \
'length($5)>1 || $6=="-" {$5="I"; $6="D"} $5=="-" || length($6)>1 {$5="D"; $6="I"} {print $1,$2,$3,$4,$5,$6}' > chr$i/$set.gmap.bim
# clean up just in case
/bin/rm chr$i/$set.sample chr$i/$set.haps
# phase haplotypes using hapi-ur
$hpr/hapi-ur -g chr$i/$set.bed -s chr$i/$set.gmap.bim -i chr$i/$set.fam -w 73 --impute -o chr$i/$set
# convert hapi-ur output for use with beagle
(echo "##fileformat=VCFv4.1";
echo "##FORMAT=<id enotype="" umber="1,Type=String,Description=\">";
echo -en "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t";
awk 'NR>2 {printf $1"."$2"\t"}' chr$i/$set.sample | sed 's/\t$/\n/g';
awk -v c=$i '{printf c"\t"$3"\t"$2"\t"$4"\t"$5"\t.\tPASS\t.\tGT";
for (i=6; i<NF; i+=2) printf "\t"$i"|"$(i+1); printf "\n"}' chr$i/$set.haps) |
java -Xmx3g -jar $bgl/b4.r1185.jar ibd=true gt=/dev/stdin usephase=true burnin-its=0 phase-its=0 out=chr$i/$set
# convert beagle output for use with DASH
awk '{sub("\\\."," ",$1); sub("\\\."," ",$3); print $1"."$2-1,$3"."$4-1,$6,$7}' chr$i/$set.hbd chr$i/$set.ibd |
$dsh/dash_cc chr$i/$set.fam chr$i/$set
# create map file with IBD clusters
cut -f1-5 chr$i/$set.clst | awk -v c=$i '{ print c,c"_"$1,$2,$3,$4,$5 }' > chr$i/$set.cmap
# create map file with IBD clusters for use with plink
awk -v c=$i '{ print c,$2,0,int(($3+$4)/2) }' chr$i/$set.cmap > chr$i/$set.map
# convert output to plink binary format
p-link --maf 0.001 --file chr$i/$set --make-bed --recode --out chr$i/$set
done
</id></pre>
<br />
You can then join the results in a single plink binary file using the code discussed in a <a href="http://apol1.blogspot.com/2013/10/merging-plink-binary-files.html">previous post</a>:
<br />
<pre># join separate autosome outputs into a single plink binary file
for i in {1..22}; do cat chr$i/$set.cmap; done > $set.dash.cmap
cp $set.fam $set.dash.fam
for i in {1..22}; do cat chr$i/$set.bim; done > $set.dash.bim
(echo -en "\x6C\x1B\x01"; for i in {1..22}; do tail -c +4 chr$i/$set.bed; done) > $set.dash.bed
</pre>
Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-2735750853096182491.post-13172011346326562302013-10-30T14:11:00.001-07:002013-10-30T14:11:53.683-07:00Longest homopolymer in coding sequence in the human genomeIt is expected that most homopolymer sequences be negatively selected in coding sequence, due to the higher mutation rate. So what is the longest homopolymer in coding sequence in the human genome? Here a short script to answer this question.<br />
<br />
First you will need a few programs:
<br />
<pre>sudo apt-get install wget samtools bedtools
</pre>
<br />
And you will need a copy of the hg19 genome:
<br />
<pre>wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/chromFa.tar.gz
tar xfO chromFa.tar.gz > hg19.fa
samtools faidx hg19.fa
</pre>
<br />
Finally, generate a bed interval file with the coordinates of the coding exons and look for homopolymers within these:
<br />
<pre><length a="" from="$7" i=""><span b="" i="" if="" to="">wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/refGene.txt.gz
zcat refGene.txt.gz | awk '{split($10,a,","); split($11,b,","); for (i=1; i<length(a); i++) {from=$7>a[i]?$7:a[i]; to=$8<b[i]?$8:b[i]; if (to>from) print $3"\t"from"\t"to}}' | bedtools sort | bedtools merge > refGene.bed
bedtools getfasta -fi hg19.fa -bed refGene.bed -fo refGene.fa
for bp in A C G T; do hp=$(for i in {1..13}; do echo -n $bp; done); grep -B1 $hp refGene.fa; done<b>
</b></span></length></pre>
<br />
It turns out that the longest homopolymer in coding sequence is 13bp long and is within the CAMKK2 gene.Unknownnoreply@blogger.com0tag:blogger.com,1999:blog-2735750853096182491.post-43571538524038856662013-10-11T14:38:00.002-07:002013-10-11T14:39:04.269-07:00Merging plink binary filesIf you have a bunch of plink binary files for the same set of people split over different chromosomes, how do you merge them? A simple answer would be to use the --merge-list option in plink. However, this will require some setting up to do and it turns out there is a much easier way to perform this.<br />
<br />
Assuming that the order of the samples within the plink files is the same (i.e. the fam files are all identical) and that each plink binary file is in SNP-major format (this is important! see here: http://pngu.mgh.harvard.edu/~purcell/plink/binary.shtml) then the following code will just work:<br />
<br />
<pre>cat plink.chr1.fam > plink.fam
for chr in {1..22} X Y; do cat plink.chr$chr.bim; done > plink.bim
(echo -en "\x6C\x1B\x01"; for chr in {1..22} X Y; do tail -c +4 plink.chr$chr.bed; done) > plink.bed</pre>
<br />
This assumes your plink files you want to merge have names like plink.chr1.bed, plink.chr2.bed, ..., plink.chrY.bed and each file corresponds to a different chromosome for your dataset. This is a typical scenario especially if you are converting a big sequencing project delivered in multiple VCF files into plink.<br />
<br />
A plink binary file in SNP-major format is just a big table preceded by three bytes that identify the file type and the order (SNP-major or individual-major). By default plink will generate SNP-major format binary files, which are ordered by SNPs, like VCF files, so you can just concatenate them after removing the three bytes headers.Unknownnoreply@blogger.com3tag:blogger.com,1999:blog-2735750853096182491.post-81529332958691518022013-10-08T10:09:00.000-07:002013-12-12T14:32:03.560-08:00Convert VCF files to PLINK binary formatHow many times have you needed to convert a VCF file to PLINK binary format? The 1000 Genomes project has a recommended tool (http://www.1000genomes.org/vcf-ped-converter) but it only works when converting small regions. Another easy to use tool is SnpSift (http://snpeff.sourceforge.net/SnpSift.html#vcf2tped).
<br />
<br />
Overall, it really isn't that complicated to convert a VCF file. The idea is to first convert it to tped format and then let plink do the job of converting that to binary format.
<br />
<br />
<pre>sudo apt-get install plink tabix
wget http://broadinstitute.org/~giulio/vcf2plink/vcf2tfam.sh
wget http://broadinstitute.org/~giulio/vcf2plink/vcf2tped.sh
chmod a+x vcf2tfam.sh vcf2tped.sh</pre>
<br />
The tools are fairly simple and they are supposed to be flexible. Now, let's suppose your VCF file is bgzip-comrpessed. A good idea might be to split it in chromosomes and generate one plink file for each chromosome. Here some code that will do that:
<br />
<br />
<pre>for chr in {1..22} X Y MT; do
tabix -H $vcf | ./vcf2tfam.sh /dev/stdin > gpc.chr$chr.tfam
tabix $vcf $chr | ./vcf2tped.sh /dev/stdin > gpc.chr$chr.tped
plink --tfile gpc.chr$chr --make-bed --out gpc.chr$chr
/bin/rm gpc.chr$chr.tfam gpc.chr$chr.tped gpc.chr$chr.nosex gpc.chr$chr.log
done</pre>
Unknownnoreply@blogger.com1tag:blogger.com,1999:blog-2735750853096182491.post-59488356176510736182013-08-28T14:39:00.000-07:002016-08-09T06:34:17.859-07:00Impute APOE and APOL1 with 23andMeIf you have paid to get your genome genotyped by 23andMe, you probably know that there are a lot of variants that are not directly genotyped by 23andMe. To obviate this problem, large datasets like the 1000 Genomes Project can be used to compare your genotypes with the genotypes of other individuals and guess what your missing genotypes are. This process is called imputation and is routinely practiced by many researchers in my area. Although, little is available online with regard to imputation of your own 23andMe data. The best link I could find was <a href="http://www.genomesunzipped.org/2013/03/learning-more-from-your-23andme-results-with-imputation.php">here</a> but the instructions provided to perform the imputation still require a fair amount of knowledge on how to use some specific tools.<br />
<br />
I think this is unacceptable and I decided to come up with my own scripts and provide a couple of examples to show how to impute your own rs429358 APOE genotype (the famous Alzheimer variant, which you will have genotyped only if you bought the more recent 23andMe v3 chip) and also how to impute the rs73885319 rs60910145 rs71785313 APOL1 genotypes (those associated to non-diabetic kidney disease in people of African descent, as shown in my previous <a href="http://dx.doi.org/10.1126/science.1193032">paper</a>) which unfortunately are still not tested by 23andMe.<br />
<br />
I wrote a few simple scripts that will allow you to impute these genotypes in a matter of minutes (or seconds, depending on how fast you can cut and paste the code below). The following instructions should work on any recent Linux Ubuntu box. You should be able to use this code on a Mac as well, provided you have the necessary basic tools already installed.<br />
<br />
The first preliminary step is to download your raw 23andMe data. Go to your <a href="https://www.23andme.com/you/download/">account</a> and download the "All DNA" data set.<br />
<br />
The second preliminary step is to install a few UNIX-friendly programs to manipulate genetic datasets, the imputation software <a href="http://faculty.washington.edu/browning/beagle/b4.html">Beagle</a>, a couple of scripts I wrote, and a template VCF file for the variants tested by 23andMe. The following commands will do:<br />
<br />
<pre>sudo apt-get install bedtools dos2unix openjdk-8-jre tabix unzip wget
wget https://faculty.washington.edu/browning/beagle/beagle.16Jun16.7e4.jar
wget http://personal.broadinstitute.org/giulio/23andme/make_vcf_template.sh
wget http://personal.broadinstitute.org/giulio/23andme/23andme_to_vcf.sh
wget http://personal.broadinstitute.org/giulio/23andme/vcf_impute.sh
wget http://personal.broadinstitute.org/giulio/23andme/23andme.vcf.gz
chmod a+x make_vcf_template.sh 23andme_to_vcf.sh vcf_impute.sh</pre>
<br />
The 23andme.vcf.gz file is a VCF template I pre-built for 23andMe v2 and 23andMe v3 chips. If you want to build it yourself (optional, you don't need to do this!), from a directory where you have stored your 23andMe raw files, run:<br />
<br />
<pre>wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b147_GRCh37p13/VCF/00-All.vcf.gz
./make_vcf_template.sh 00-All.vcf.gz genome_*_Full_*.zip</pre>
<br />
This step takes time as the whole dbSNP dataset (~1.4GB) will need to be downloaded and processed.<br />
<br />
Now the first step is to convert the raw downloaded data to a more standard format, the <a href="http://www.1000genomes.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-41">VCF</a> format:<br />
<br />
<pre>./23andme_to_vcf.sh 23andme.vcf.gz genome_*_Full_*.zip</pre>
<br />
Where genome_*_Full_*.zip should be substituted with the name of the zip file you downloaded from 23andMe. The second and last step is to create a list of variants (from a single locus) to be imputed and pass the genomic locations of these variants to a script which will run the imputation.<br />
<br />
The script provided will automatically download a piece of the 1000 Genomes Project reference panel around the variants of interest, subset this genotype panel to the variants to be imputed and the variants available in your VCF file, and run the imputation with this minimal reference panel. This way the imputation process will be immediate and you will have your results right away.<br />
<br />
For APOE, we are only interested in the rs429358 variant and the rs7412 variant (which is already genotyped both in the 23andMe v2 and the 23andMe v3 chip), so the code will look like this:<br />
<br />
<pre>echo -e "19\t45411940\t45411941\n19\t45412078\t45412079" > apoe.bed
./vcf_impute.sh beagle.16Jun16.7e4.jar apoe.bed genome_*_Full_*.vcf.gz</pre>
<br />
Where genome_*_Full_*.vcf.gz should be substituted with the name of the compressed VCF file generated in the previous step. For APOL1 we are interested instead in three variants:<br />
<br />
<pre>echo -e "22\t36661905\t36661906\n22\t36662033\t36662034\n22\t36662040\t36662047" > apol1.bed
./vcf_impute.sh beagle.16Jun16.7e4.jar apol1.bed genome_*_Full_*.vcf.gz</pre>
<br />
The results will be displayed on screen in VCF file format with best guess genotypes and individual genotype likelihoods. Of course, the results are only probabilistic, so you will need to interpret them accordingly. For the biological interpretation instead, snpedia has a good summary, both for <a href="http://snpedia.com/index.php/APOE">APOE</a> and for <a href="http://snpedia.com/index.php/APOL1">APOL1</a>.<br />
<br />
If you have any comments, complaint, or suggestions, please send me an <a href="mailto:giulio.genovese@gmail.com">email</a>! No perl code was used for any of the scripts.Unknownnoreply@blogger.com7