基因组survey调研图分析

  • 杂合度：大多数真核生物有两套匹配的染色体组（同源染色体） ; 也就是说，它们是二倍体。二倍体生物，如果两个同源染色体上相同基因座上两个等位基因相同，则该生物在该位点是纯合的。如果不同，则该生物体在该基因座处是杂合的。

  • 一个基因的两个不同等位基因（一个野生型等位基因和一个突变等位基因）时，它在基因座处是杂合的。细胞或有机体被称为杂合子， 专门用于所讨论的等位基因，因此，杂合性是指特定的基因型。杂合基因型由大写字母（代表显性/野生型等位基因）和小写字母（代表隐性/突变等位基因）表示，例如“ Rr”或“ Ss”。大写字母通常是先写的。

  • 杂合基因型可以比纯合显性基因型或纯合隐性基因型具有更高的相对适应度-这被称为杂合子优势。

  • 降低杂合度可以通过很多代近交来实现。

一、前言

对于高等真核生物（特别是高等植物）来讲，它们的基因组一般比较复杂，通常伴随着较高的杂合度或者重复序列等，进行基因组De novo也并非一件容易的事情。特别是进行高精细基因组图谱的构建时，只基于常规DNA文库的二代测序根本无法满足需要，还需结合三代测序、BioNano、Hi-C数据等才可获得高质量的基因组。相对而言，这时的测序成本也会是相当大的一笔开销，需要谨慎对待。因此在对高等真核生物进行全基因组De novo之前，我们需要设法提前获知该物种基因组特征的一些信息，为后续的测序方案、基因组组装方案、基因组结构注释等提供参考依据。这种情况下，我们一般会在基因组大规模测序或者正式组装之前，首先构建DNA小片段文库进行中低深度的二代测序，使用PE文库测序所得的reads信息进行基因组Survey分析以初步评估基因组特征，这是很重要的一步。

因此，在对一个新物种测序时，首先需要对该物种的基因组结构有所了解，包括以下几个内容：

（1）基因组大小
　　基因组大小的获取，关系到对以后组装结果大小的准确判断与否；基因组太大（>10Gb），可能会超出了目前denovo组装基因组软件的对机器存储的要求，从客观条件上讲是无法实现组装的。一般物种的基因组大小可以从公共数据库查到。如果没有搜录，需要考虑通过实验（流式细胞仪福尔根 染色/定量pcr/）或Kmer估计法来获得基因组大小。

（2）杂合度大小
　　杂合度对基因组组装的影响主要体现在不能合并姊妹染色体。杂合度高的区域，会把两条姊妹染色单体都组装出来，从而造成组装的基因组偏大于实际的基因组大小。一般是通过SSR在测序亲本的子代中检查SSR的多态性。杂合度如果高于0.5%，则认为组装有一定难度。杂合度高于1%则很难组装出来。杂合度高，并不是说组装不出来，而是说，装出来的序列不适用于后续的生物学分析。比如拷贝数、基因完整结构。

（3）重复片段大小

总结：基因组越大，杂合度也大，重复片段越大，该物种的组装难度就越大。通常我们会通过genome survery分析，对以上几个指标进行简单评估，核心就是通过kme 分布来进行评估。

二、Kmer

K-mer是对于我们拿到的测序数据迭代选取长度为K的序列片段，示意图如下：

K-mer原理图

图中第一行为提供的序列，下方为K=5时取出的K-mer序列。

拿到所有的k-mer序列后，就可以根据k-mer序列出现的频率进行基因组复杂度评估了。

三、Kmer分析

推荐使用Jellyfish+GenomeScope的组合进行基因组K-mer分析。

1. Jellyfish

安装jellyfish

conda install -c bioconda jellyfish

使用jellyfish软件产生的kmer频数分布数据

jellyfish count  -m 21 -s 20G -t 20 -o 21mer_out  -C  <(zcat test_1.fq.gz) <(zcat test_2.fq.gz)
jellyfish histo -t -o 21mer_out.histo 21mer_out

Options:
　　-m kmer length；
　　-s memory (-s) , Hash大小, 根据文件大小确定；
　　-t threads 线程；
　　-o 输出前缀；
　　-C 统计正负链，一定要用；

注意事项：

  1. 绝对不要用--min-qual-char或其他参数，它们会将低质量的碱基替换成N
  2. 在测序时由于不知道测得到底是DNA的哪一条链，因此k-mer及其互补链其实是等价的，所以一定要用-C参数

-C，canonical是什么含义？

  • 对于 fasta 类型数据来说，这只是一条单链的DNA序列，本身就不存在的反向重复，因此在对其进行分析的时候，没有必要使用-C参数。
  • 对于 fastq 类型数据来说，由于测序时将基因组DNA（当然也可以是RNA）打断之后，对2条链都进行了测序。但是在进行基因组大小估计、选择亚基因组特异性Kmer等分析时，都是默认一条DNA序列的情况，因此在针对fastq类型的数据进行分析时，我们就 需要加上-C参数 。

多个lane时，可以＊匹配fq文件

$ jellyfish count -C -m 21 -s 1000000000 -t 10 *.fastq -o 21mer_out

将数据导入R语言中，进行作图

pdf("21_mer.out.pdf")
dataframe <- read.table("21mer_out.histo")
plot(dataframe[1:200,], type="l")
dev.off()
image

由于只有一个主峰，说明该物种的杂合度并不高，基本上也就是二倍体。如果图中出现多个峰，说明它可能是多倍体或者是基因组杂合度高。

基因组大小(G)估计算法为:


其中Kdepth为K-mer的期望测序深度,Knum为K-mer的总数。 通常将K-mer深度分布曲线的峰值作为其期望深度。

2. GenomeScope根据Jellyfish统计出来的K-mer分布结果进行统计预估

将Jellyfish得到的reads.histo 传给genomescope等待几分钟后，就会得到结果

基因组的杂合性和使得来自杂合片段的K-mer深度较纯合区段降低50%。如果目标基因组有一定的杂合性，会在k-mer深度分布曲线主峰位置(c)的1/2处(c/2)出现一个小峰。杂合度越高，该峰越明显。

注意：
基因组预测大小和第一个页面中 Max kmer coverage 密切相关。GenomeScope默认会过滤掉出现1000次以上的kmers，避免细胞器基因组的影响，如果你觉得基因组小了，那么就把数值调整的大一点。

四、其他

1. 基于组装

基于K-mers可以较好的预测基因组大小，并定性的了解基因组的复杂情况，如果想更具体的了解基因组的复杂度，可以先将50X以上的段片段进行组装，然后进行分析。

组装的工具比较多，推荐用SOAPdenovo，因为速度快。

新建一个contig.config, 增加如下内容

max_rd_len=150
[LIB]
avg_ins=200
reverse_seq=0
asm_flags=3
rd_len_cutoff=100
rank=1
pair_num_cutoff=3
map_len=32
q1=read_1.fq
q2=read_2.fq

组装出参考序列

~/opt/biosoft/SOAPdenovo2/SOAPdenovo-63mer all -s contig.config -R -K 63 -p 30 -o assembly/graph

最后graph.scafSeq是拼接后的序列, 提取出大于300bp的序列.

# adjust format
bioawk -c fastx -v name=1 '{if(length($seq)>300) print ">"name "\n" $seq;name+=1}' assembly/graph.scafSeq >contig.fa

2.杂合度估计

将原来的序列回贴到contig上，并用samtools+bcftools进行snp calling.统计变异的碱基占总体的比例。

mkdir -p index
bwa index contig.fa -p index/contig
bwa mem -v 2 -t 10 index/contig read_1.fq read_2.fq | samtools sort -n > align.bam
samtools mpileup -f contig align.bam | bcftools call -mv -Oz -o variants.gz

一方面由于SOAPdenovo组装过程中会出错, 另一方面samtools在变异检测上也存在很高的假阳性, 所以总得先按照深度和质量过滤一批假阳性。

bcftools view -i ' DP > 30 && MQ > 30' -H variants.vcf.gz | wc -l
# 325219, 无过滤是445113

变异数目占基因组大小的比例就是杂合度。我的contig大概是200M，找到0.3M左右的变异，也就是0.0015，即0.15%.

3. 重复序列估计

基于同源注释，用RepeatMasker寻找重复序列. 这里要注意分析的fasta的ID不能过长，也就是最好是>scaffold_1这种形式，不然会报错。

~/opt/biosoft/RepeatMsker/RepeatMasker -e ncbi -species arabidopsis -pa 10 -gff -dir ./ contig.fa
# -e ncbi
# -species 选择物种 用~/opt/biosoft/RepeatMasker/util/queryRepeatDatabase.pl -tree 了解
# -pa 并行计算
# -gff 输出gff注释
# -dir 输出路径

输出结果中主要关注如下三个

  • output.fa.masked, 将重复序列用N代替
  • output.fa.out.gff, 以gff2形式存放重复序列出现的位置
  • output.fa.tbl, 该文件记录着分类信息
==================================================
file name: anno.fasta
sequences:         62027
total length:  273135210 bp  (273135210 bp excl N/X-runs)
GC level:         36.80 %
bases masked:   79642191 bp ( 29.16 %)
==================================================

也就是说我们的物种有30%的重复序列，作为参考，拟南芥125Mb 14%重复序列, 水稻389M，36%重复

软件安装
安装RepeatMasker

cd ~/src
wget http://tandem.bu.edu/trf/downloadstrf409.linux64
mv trf409.linux64 ~/opt/bin/trf
chmod a+x ~/opt/bin/trf
# RMBlast
cd ~/src
wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.6.0/ncbi-blast-2.6.0+-src.tar.gz
wget http://www.repeatmasker.org/isb-2.6.0+-changes-vers2.patch.gz
tar xf ncbi-blast-2.6.0+-src
gunzip isb-2.6.0+-changes-vers2.patch.gz
cd ncbi-blast-2.6.0+-src
patch -p1 < ../isb-2.6.0+-changes-vers2.patch
cd c++
./configure --with-mt --prefix=~/opt/biosoft/rmblast --without-debug && make && make install
# RepeatMasker
cd ~/src
wget http://repeatmasker.org/RepeatMasker-open-4-0-7.tar.gz
tar xf RepeatMasker-open-4-0-7.tar.gz
mv RepeatMasker ~/opt/biosoft/
cd ~/opt/biosoft/RepeatMasker
## 解压repbase数据到Libraries下
## 配置RepatMasker
perl ./configure
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 157,298评论 4赞 360
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 66,701评论 1赞 290
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 107,078评论 0赞 237
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 43,687评论 0赞 202
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 52,018评论 3赞 286
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 40,410评论 1赞 211
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 31,729评论 2赞 310
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 30,412评论 0赞 194
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 34,124评论 1赞 239
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 30,379评论 2赞 242
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,903评论 1赞 257
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 28,268评论 2赞 251
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 32,894评论 3赞 233
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 26,014评论 0赞 8
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,770评论 0赞 192
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 35,435评论 2赞 269
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 35,312评论 2赞 260