【直播】我的基因组52:X和Y染色体的同源区域探索

很久以前,我其实就遇到过通过NGS测序数据来判定性别的难题(搜索我博客即可查看详情),本次探究自己的基因组得到的统计结果与常识不符,所以我可以肯定是我们的常识太浅显了。

【直播】我的基因组48:我可能测了一个假的全基因组

【直播】我的基因组49:Y染色体的SNV不能用常规流程来找?

【直播】我的基因组50:从测序深度和位点间距来看SNV分布情况

通过自己的测序数据的详细分析,我才知道PAR(pseudoautosomal region)。这样的X,Y染色体大量同源,说到底是测序片段压根无法准确定位,所以说所谓的X,Y染色体是单倍体的常识,在这里完全错误的。这些区域目前有29个基因,那么对这29个基因来说,其实就跟定位在常染色体上一样,有两个拷贝的!

这些区域在hg38的参考基因组坐标如下;

The locations of the PARs within GRCh38 are:

PAR1: chrY:10,000-2,781,479 and chrX:10,000-2,781,479 [7]

PAR2: chrY:56,887,902-57,217,415 and chrX:155,701,382-156,030,895 [7]

PAR3: chrY:3,571,959-5,881,959 and chrX:89,145,000-92,745,001 [3]

那么我们就可以通过自己的数据处理能力来探索一下X和Y染色体的同源区有多少,是哪里的问题!

首先下载X,Y染色体的fasta序列,在UCSC上面下载即可。

然后把X染色体构建bwa的索引。

接着模拟一个Y染色体的测序数据,模拟的程序很简单,模拟Y染色体的测序片段(PE100,insert400)。

最后把模拟测序数据比对到X染色体的参考,统计一下比对结果即可!

我自己看sam文件也发现真的同源性好高呀,总共就模拟了380万reads,就有120万是百分百比对上了。

所以对女性个体来说,测序判断比对到Y染色体是再正常不过的了。如果要判断性别,必须要找那些X,Y差异性区段!对男性来说,更是如此!

本次测试涉及到的文件如下:

shell脚本如下:

  1. cd tmp/chrX_Y/hg19/
  2. wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrX.fa.gz;
  3. wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chrY.fa.gz;
  4. gunzip chrX.fa.gz
  5. gunzip chrY.fa.gz
  6. ~/biosoft/bwa/bwa-0.7.15/bwa index chrX.fa
  7. ~/biosoft/bwa/bwa-0.7.15/bwa mem -t 5 -M chrX.fa read*.fa >read.sam
  8. samtools view -bS read.sam >read.bam
  9. samtools flagstat read.bam
  10. samtools sort -@ 5 -o read.sorted.bam read.bam
  11. samtools view -h -F4 -q 5 read.sorted.bam |samtools view -bS|samtools rmdup - read.filter.rmdup.bam
  12. samtools index read.filter.rmdup.bam
  13. samtools mpileup -ugf ~/tmp/chrX_Y/hg19/chrX.fa read.filter.rmdup.bam |bcftools call -vmO z -o read.bcftools.vcf.gz

对Y染色体随机抽取模拟测序片段的程序如下(这个程序我不想给文字版的,希望大家可以自己手动敲一遍,在我们的生信技能树论坛上面提交自己的感悟:http://www.biotrainee.com/thread-696-1-1.html):

这个测序待改进的地方太多了,比如可以过滤掉N含量过多的片段(我只是把全部是N的地方去除了),可以设置插入片段为参数,而且打断的片段不应该是稳定的600bp,而且可以改成PE150的测序,或者更长,模拟一下看看是不是3代测序的超长片段,就能解决这个问题。

建bwa索引的log日志如下:

仔细打开比对结果sam文件可以继续探索,有不少比对结果含义XA:Z,说明即使是这100个碱基在X染色体也有多个定位!

【直播】我的基因组(十三):了解sam格式比对结果

甚至对这个sam文件可以做variation的calling,然后放到IGV里面去看看!

最后找到的variation也可以统计一下:

96180个 0/1

181020 个1/1

当然,这里我模拟的是4X 的数据,所以找到的variation不会太准确,但是我模拟的精确数据,其实不应该有杂合的variation,但结果还是有一些~

毕竟这种比对也太诡异了,看来我对BWA软件的理解还不够透彻!

请参与本次直播基因的同学继续我的思路探索下去,模拟PE150,甚至miseq的PE250的测序片段看看比对情况如何,或者模拟三代测序仪的。

还可以下载hg38参考基因组的X,Y序列,只有你实践的越多你才能学到更多!

只有你实践的越多你才能学到更多!

只有你实践的越多你才能学到更多!

只有你实践的越多你才能学到更多!

参考:https://en.wikipedia.org/wiki/Pseudoautosomal_region

文:Jimmy

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏听雨堂

振幅和成交量的关系

用广晟有色的历史数据,用sklearn进行回归,数据如下: ? 假设每日振幅和成交量以及价格是有关系的,于是构造: # coding=utf-8 from pa...

2008
来自专栏生信技能树

【直播】我的基因组47:测序深度和GC含量的关系

在前面我们用 ChIP-seq 的分析方法可视化了一下我的 WGS数据,结果我们的测序深度分布居然是跟基因组的genomic feature相关。 比如在TSS...

38111
来自专栏生信宝典

生信宝典之傻瓜式(四)蛋白蛋白互作网络在线搜索

傻瓜系列重启了,今天要介绍的是一款在线查询蛋白-蛋白互作网络的工具 STRING (https://string-db.org/)。 STRING数据库收录了2...

3735
来自专栏吉浦迅科技

DAY81:阅读Compute Capability 5.x

When a multiprocessor is given warps to execute, it first distributes them among...

1383
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 3

前面两篇文章(高通量数据分析必备|基因组浏览器使用介绍 - 1和高通量数据分析必备|基因组浏览器使用介绍 - 2)介绍了EPGG的基本使用、各部分特征、Trac...

1115
来自专栏生信技能树

【直播】我的基因组81:看看我的vcf文件的vaf分布情况

这一讲中,我们对VCF中的"VAF"简单的来看一起,如果你对VCF文件还不了解的话,那你就要自我批评一下了。在基因组直播刚开始的时候,我还专门对VCF文件进行了...

5716
来自专栏数据小魔方

R语言可视化——图表美化与套用主题(下)

昨天的分享跟大家简单介绍了关于柱形图图表元素美化的思路,今天接着分享关于套用主题。 因为单独使用代码来调整单个图表元素,实在是太费劲了,更何况图表的细节元素有那...

2766
来自专栏Y大宽

从原始芯片.cel数据到权重基因共表达网络(WGCNA)详细流程

看这个之前,可以先看WGCNA的一些理论背景知识 看完整个之后可以去看WGCNA关键模块和hub基因筛选

3803
来自专栏生信宝典

高通量数据分析必备|基因组浏览器使用介绍 - 1

基因组浏览器是高通量测序分析的一个重要的可视化工具。相比于最终提供的表格,基因组浏览器可以提供更多的信息,如直观展示突变位点、查看有无新转录本或新的可变剪接形式...

1222
来自专栏生信技能树

【直播】我的基因组59:CNV初步探索

好久不见,基因组直播又来了。这篇推送是对SNV进行一个初步探索。 单纯的一个样本来找CNV,总是不太准确的,但还是那句话,毕竟是自己的基因组,硬着头皮也要上。当...

39613

扫码关注云+社区