Phylogenomic_Tutorial || Ortholog Detection

Tutorial

[TOC]

我们在做物种尺度上的系统发育推断时候，通常是利用基因树gene tree来代表物种树species tree。而这个推断的大的前提是基因的分化历史与物种的分化历史完全相同。而实际上一些基因会由于某物种内自身的复制、丢失，谱系分选等原因使得基因树不能完全代表物种树。这时候一种方法就是利用多个同源基因orthologs联合建树的方法去做物种水平上的系统推断。本篇教程即是为了从不同物种或转录组中获得一对一的直系同源基因。

Preparation

软件：

  • BLAST
  • PAML 利用其中的codeml计算 dn/ds速率。
  • MAFFT 多序列比对
  • BMGE 保守区段的筛选
  • find_ortholog.py十分实用的脚本。

数据：

  • 利用预先的2231个exon的query序列从12个有genome或transcriptome的物种筛选直系同源基因。query序列大致是属于没有复制/缺失的序列。12个物种是包括11个属于正在快速分化的cichlid鱼，和外类群Zebrafish
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
## 下载10个物种的基因组

# Download the genome assembly of Amphilophus citrinellus.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/751/415/GCA_000751415.1_Midas_v5/GCA_000751415.1_Midas_v5_genomic.fna.gz
  gunzip GCA_000751415.1_Midas_v5_genomic.fna.gz
  mv GCA_000751415.1_Midas_v5_genomic.fna ampcit.fasta
  
  # Download the genome assembly of Andinoacara coeruleopunctatus.
  wget http://cichlid.gurdon.cam.ac.uk/Andinoacara_coeruleopunctatus_final_min1000bp_scaffolds.fa.gz
  gunzip Andinoacara_coeruleopunctatus_final_min1000bp_scaffolds.fa.gz
  mv Andinoacara_coeruleopunctatus_final_min1000bp_scaffolds.fa andcoe.fasta
  
  # Download the genome assembly of Oreochromis niloticus.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/vertebrate_other/Oreochromis_niloticus/all_assembly_versions/GCF_001858045.1_ASM185804v2/GCF_001858045.1_ASM185804v2_genomic.fna.gz
  gunzip GCF_001858045.1_ASM185804v2_genomic.fna.gz
  mv GCF_001858045.1_ASM185804v2_genomic.fna orenil.fasta
  
  # Download the genome assembly of Astatotilapia burtoni.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/239/415/GCF_000239415.1_AstBur1.0//GCF_000239415.1_AstBur1.0_genomic.fna.gz
  gunzip GCF_000239415.1_AstBur1.0_genomic.fna.gz
  mv GCF_000239415.1_AstBur1.0_genomic.fna astbur.fasta
  
  # Download the genome assembly of Metriaclima zebra.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/238/955/GCF_000238955.4_M_zebra_UMD2a//GCF_000238955.4_M_zebra_UMD2a_genomic.fna.gz
  gunzip GCF_000238955.4_M_zebra_UMD2a_genomic.fna.gz
  mv GCF_000238955.4_M_zebra_UMD2a_genomic.fna metzeb.fasta
  
  # Download the genome assembly of Pundamilia nyererei.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/239/375/GCF_000239375.1_PunNye1.0//GCF_000239375.1_PunNye1.0_genomic.fna.gz
  gunzip GCF_000239375.1_PunNye1.0_genomic.fna.gz
  mv GCF_000239375.1_PunNye1.0_genomic.fna punnye.fasta
  
  # Download the genome assembly of Neolamprologus brichardi.
  wget ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/239/395/GCF_000239395.1_NeoBri1.0//GCF_000239395.1_NeoBri1.0_genomic.fna.gz
  gunzip GCF_000239395.1_NeoBri1.0_genomic.fna.gz
  mv GCF_000239395.1_NeoBri1.0_genomic.fna neobri.fasta
  
  # Download the genome assembly of Neolamprologus marunguensis.
  wget http://evoinformatics.eu/neo_assemblies/Ma.scf.fasta.gz
  gunzip Ma.scf.fasta.gz
  mv Ma.scf.fasta neomar.fasta

  # Download the genome assembly of Neolamprologus gracilis.
  wget http://evoinformatics.eu/neo_assemblies/Gr.scf.fasta.gz
  gunzip Gr.scf.fasta.gz
  mv Gr.scf.fasta neogra.fasta

  # Download the genome assembly of Neolamprologus olivaceous.
  wget http://evoinformatics.eu/neo_assemblies/Ol.scf.fasta.gz
  gunzip Ol.scf.fasta.gz
  mv Ol.scf.fasta neooli.fasta

### Outgroup reference外类群
  # Download the genome assembly of Danio rerio.
  wget http://ftp.ensembl.org/pub/release-87/fasta/danio_rerio/dna/Danio_rerio.GRCz10.dna.toplevel.fa.gz
  gunzip Danio_rerio.GRCz10.dna.toplevel.fa.gz
  mv Danio_rerio.GRCz10.dna.toplevel.fa danrer.fasta

对于没有基因组的物种Ophthalmotilapia ventralis，可以在NCBI的Nucleotide数据库找mRNA数据。https://www.ncbi.nlm.nih.gov/nucleotide/

1
mv sequence.fasta ophven.fasta

直系同源基因的初筛Identifying putative orthologs

用作query的exon.fasta 文件中的序列为蛋白pep序列。因为在用于相对较大尺度的物种间时，其相较于nucl序列能更好的去鉴别多的orthologs。

  • 构建12个genome/transcriptome的blast database
1
2
3
4
for i in ??????.fasta
  do
        makeblastdb -in ${i} -dbtype nucl
done

单个序列query寻找所有物种中的orthologs

  1. 我们首先做个测试利用query中的一个序列去操作，利用find_prtholog.py

    1
    2
    3
    4
    5
    6
    # 1. 测试数据ENSDARE00000125484序列
    head -n 2 exons.fasta > tmp.fasta

    # 2. 根据find_prtholog.py脚本
    ### -t 表示query为pep序列，对subject使用tblastn程序
    python3 find_orthologs.py -t tmp.fasta orenil.fasta

  2. 批量对12个genome进行操作。

1
2
3
4
5
6
7
8
#1. 将已经makeblastdb的基因组fasta文件输出到文件
ls danrer.fasta > subjects.txt
ls ??????.fasta | grep -v danrer >> subjects.txt

#2. 对query序列比对到多个基因组的subject 基因组上找orthologs
### -s 1 对应 BLAST / TBLASTN option -culling_limit 1参数；比对的严格程度
### -r/--refine 表示找到的序列进行mafft多序列比对。
python3 find_orthologs.py -t -s 1 --refine tmp.fasta subjects.txt

nucl.fasta文件会进行mafft多序列比对，其中ophven没找见也不代表没有，可能因为bitscore值过低而被findortholog.py程序所删除

多个query序列寻找所有物种中的orthologs

  1. 利用filter_queries.rb对query序列进行初筛
1
2
3
4
### filter_queries.rb 输入.fasta 输出.fasta 过滤序列长度小于60 过滤bitscore低于51
ruby filter_queries.rb exons.fasta exons_red.fasta 60 51

### 结果显示：Removed 1696 out of 2231 sequences.
  1. 对此1696的query seqs进行比对
1
2
3
4
5
6
7
8
9
10
11
12
#1. 利用split_queries.rb程序对1696序列分割为5个
ruby split_queries.rb exons_red.fasta 5

#2. 并行计算

for i in exons_red.??.fasta
  do
        python3 find_orthologs.py -t -s 1 --refine ${i} subjects.txt &
done
  
## 多种方法都可以实现并行运算。程序根据blast时间长短，耗时较久。
## 结果会出现1696个fasta文件，和mafft多序列比对后的文件。

过滤orthologs

blast比对并非完全能找到不同物种中的直系同源基因，还存在可能会找到不同物种中的旁系同源基因paralogs，或者比对错误的基因序列misaligned。我们需要将这些可能错误的序列进行过滤，从而能更好的用于phylogeny以及分化时间推断。

过滤1：根据bitscore值进行过滤

因为我们的query是用相对较远的zebrafish进行blast，所以在cichlid类群中找到的orthologs的bitscore值应该是非常相近的，而如果并非如此，过低的bitscore值就可能代表paralog或错误blast。因此就可以直接根据此bitscore设定阈值进行过滤。根据ruby脚本filter_sequences_by_bitscore.rb

1
2
# 包含有所有序列的输入directory, 输出的directory, 根据最大的bitscore的百分位数设置的阈值
ruby filter_sequences_by_bitscore.rb 01 02 0.9

结果会去掉文件中bitscore低于max * 0.9阈值的序列；

过滤2：根据dN/dS的选择压力阈值进行过滤

基本假设就是，如果这个cichlid类群中的序列彼此间都为orthologs，那么与外类群的zebrafish的query 序列相比都会有相似的选择信号值（以dn/ds来表示）。如果一个序列是paralog，那么此paralog可能会因为新功能化导致出现不同的选择压力值，因此会出现增高的dn/ds值。所有可以设置一个dN/dS阈值，过滤高的序列；同时，去除受选择的序列也更有利于后续分化时间的分析。

利用PAML软件中codeml来计算dn/ds值，同时也利用作者写好的filter_sequences_by_dNdS.rb 来自动化的根据阈值进行过滤。

1
2
# 包含有所有序列的输入directory, 输出的directory； dn/ds阈值0.3
ruby filter_sequences_by_dNdS.rb 02 03 0.3

结果同样会对每个文件中的omega值低于0.3的物种序列进行过滤。阈值可根据实际情况进行调整，但一个基本准则是对于系统发育分析来说，在序列/信息位点够多的情况下，remove rather too much than too little.

过滤3：利用BMGE保留信息位点

经过以上两部，我们assume物种间的序列都为orthologs，接下来再利用BMGE来筛选保守的信息位点，以便于更好的系统发育分析；由于是cds序列，需按照三联读码框的方式来过滤。利用filter_sites_with_BMGE.rb 脚本来过滤。

1
2
## 包含有所有序列的输入directory, 输出的directory；BMGE.jar程序所在位置； maximum gap rate of 0.2； maximum entropy score 0.5
ruby filter_sites_with_BMGE.rb 03 04 BMGE.jar 0.2 0.5

过滤4：去除过多缺失的序列

经过以上的过滤，会出现许多序列文件中的某一些物种是缺失的，利用程序filter_exons_by_missing_data.rb 程序4个参数。

1
2
## 输入文件夹，输出文件夹；允许的完全缺失序列数目；每条序列最短的长度
ruby filter_exons_by_missing_data.rb 04 05 1 150

表示超过1个物种都是缺失的序列舍弃；序列长度低于150bp的丢弃。

过滤5：根据GC值的标准差进行过滤

基本假设就是物种间过高或过低的GC值会导致系统发育分析的可能错误。因此我们去除这些超过sd范围的序列。这一步需要去除outgroup序列。也就是我们的query序列。根据脚本filter_exons_by_GC_content_variation.rb 进行过滤

1
2
## 过滤sd值大于0.03的序列
ruby filter_exons_by_GC_content_variation.rb 06 07 0.03

注意：即使是经过以上的过滤方法之后，还是会出现在多序列比对中，一些物种的序列存在错误的情况，因此如果序列相对较少，还是需要一条条人工的过滤那些明显错误比对的序列。

后续操作

实际上由于重组recombination或是不完全谱系分选incomplete lineage sorting，genome中的一些片段gene tree是不能代表物种树的。由于后续会进行multi-species-coalescent model下的系统推断；因此，在这里我们不会对过滤后的1000多个序列进行串联成单个super matrix。但是，一些序列由于过滤会导致过短而使得gene tree的信息位点都不够用，所以如果一些序列太短的话，一种折衷的方法是对同属于一个gene中的exon序列进行合并，或类似的小范围合并，假设这段序列之间是相同的进化历史。这就是其它的手段可以实现了。

References

Malmstrøm et al. (2016)
https://github.com/uio-cees/teleost_genomes_immune/blob/master/S11_S14_ortholog_detection_and_filtering
Species-tree inference

本文标题:Phylogenomic_Tutorial || Ortholog Detection

文章作者:Tianpeng Wang

发布时间:2021-01-21, 01:00:19

最后更新:2021-01-21, 01:05:09

原始链接:http://wangtp.top/2021/01/21/Phylogenomic_Tutorial_Ortholog%20Detection/

许可协议: "署名-非商用-相同方式共享 4.0" 转载请保留原文链接及作者。

文章目录
  1. 1. Preparation
  2. 2. 直系同源基因的初筛Identifying putative orthologs
    1. 2.1. 单个序列query寻找所有物种中的orthologs
    2. 2.2. 多个query序列寻找所有物种中的orthologs
  3. 3. 过滤orthologs
    1. 3.1. 过滤1：根据bitscore值进行过滤
    2. 3.2. 过滤2：根据dN/dS的选择压力阈值进行过滤
    3. 3.3. 过滤3：利用BMGE保留信息位点
    4. 3.4. 过滤4：去除过多缺失的序列
    5. 3.5. 过滤5：根据GC值的标准差进行过滤
    6. 3.6. 后续操作
  4. 4. References
|