【直播】我的基因组52：X和Y染色体的同源区域探索
很久以前，我其实就遇到过通过NGS测序数据来判定性别的难题(搜索我博客即可查看详情)，本次探究自己的基因组得到的统计结果与常识不符，所以我可以肯定是我们的常识太浅显了。
【直播】我的基因组49:Y染色体的SNV不能用常规流程来找？
【直播】我的基因组50:从测序深度和位点间距来看SNV分布情况
通过自己的测序数据的详细分析，我才知道PAR(pseudoautosomal region)。这样的X,Y染色体大量同源，说到底是测序片段压根无法准确定位，所以说所谓的X,Y染色体是单倍体的常识，在这里完全错误的。这些区域目前有29个基因，那么对这29个基因来说，其实就跟定位在常染色体上一样，有两个拷贝的！
这些区域在hg38的参考基因组坐标如下；
The locations of the PARs within are:
PAR1: chrY:10,000-2,781,479 and chrX:10,000-2,781,479
PAR2: chrY:56,887,902-57,217,415 and chrX:155,701,382-156,030,895
PAR3: chrY:3,571,959-5,881,959 and chrX:89,145,000-92,745,001
那么我们就可以通过自己的数据处理能力来探索一下X和Y染色体的同源区有多少，是哪里的问题！
首先下载X,Y染色体的fasta序列，在UCSC上面下载即可。
然后把X染色体构建bwa的索引。
接着模拟一个Y染色体的测序数据，模拟的程序很简单,模拟Y染色体的测序片段（PE100，insert400）。
最后把模拟测序数据比对到X染色体的参考，统计一下比对结果即可！
我自己看sam文件也发现真的同源性好高呀，总共就模拟了380万reads，就有120万是百分百比对上了。
所以对女性个体来说，测序判断比对到Y染色体是再正常不过的了。如果要判断性别，必须要找那些X,Y差异性区段！对男性来说，更是如此！
本次测试涉及到的文件如下：
shell脚本如下：
cd tmp/chrX_Y/hg19/
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz;
wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz;
gunzip chrX.fa.gz
gunzip chrY.fa.gz
~/biosoft/bwa/bwa-0.7.15/bwa index chrX.fa
~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M chrX.fa read*.fa >read.sam
samtools view -bS read.sam >read.bam
samtools flagstat read.bam
samtools sort -@ 5 -o read.sorted.bam read.bam
samtools view -h -F4 -q 5 read.sorted.bam |samtools view -bS|samtools rmdup - read.filter.rmdup.bam
samtools index read.filter.rmdup.bam
samtools mpileup -ugf ~/tmp/chrX_Y/hg19/chrX.fa read.filter.rmdup.bam |bcftools call -vmO z -o read.bcftools.vcf.gz
对Y染色体随机抽取模拟测序片段的程序如下(这个程序我不想给文字版的，希望大家可以自己手动敲一遍，在我们的生信技能树论坛上面提交自己的感悟：)：
这个测序待改进的地方太多了，比如可以过滤掉N含量过多的片段(我只是把全部是N的地方去除了)，可以设置插入片段为参数，而且打断的片段不应该是稳定的600bp，而且可以改成PE150的测序，或者更长，模拟一下看看是不是3代测序的超长片段，就能解决这个问题。
建bwa索引的log日志如下：
仔细打开比对结果sam文件可以继续探索，有不少比对结果含义XA:Z，说明即使是这100个碱基在X染色体也有多个定位！
甚至对这个sam文件可以做variation的calling，然后放到IGV里面去看看！
最后找到的variation也可以统计一下：
96180个 0/1
181020 个1/1
当然，这里我模拟的是4X 的数据，所以找到的variation不会太准确，但是我模拟的精确数据，其实不应该有杂合的variation，但结果还是有一些~
毕竟这种比对也太诡异了，看来我对BWA软件的理解还不够透彻！
请参与本次直播基因的同学继续我的思路探索下去，模拟PE150,甚至miseq的PE250的测序片段看看比对情况如何，或者模拟三代测序仪的。
还可以下载hg38参考基因组的X,Y序列，只有你实践的越多你才能学到更多！
只有你实践的越多你才能学到更多！
只有你实践的越多你才能学到更多！
只有你实践的越多你才能学到更多！
参考：
文：Jimmy
图文编辑：吃瓜群众