1. 简介
在本教程中，咱们将把两个 Heliconius 蝴蝶物种的一条染色体（包括 optix 基因）与泛基因组进行比对。
泛基因组是指一个物种的所有个别同享的完好基因组序列，以及特定个别或亚群所独有的可变基因组序列。咱们将运用 seq-seq-pan 构建泛基因组比对，运用一些自定义 Python 脚原本解析输出，并运用 R 来可视化比对。此外，咱们将把发育中的头部和翅膀组织的转座元件（TE）注释和染色质可及性图谱（ATAC-seq）的坐标转换到泛基因组坐标空间，并将它们添加到该图中。
最终成果应如下所示：
2. 输入数据
关于本练习，您能够导航至 ensembl.lepbase.org/ 并单击 Heliconiuserato demophoon (v1) 和 Heliconius melpomene melpomene (Hmel2) 基因组组件的链接。在右上角，您能够查找“optix”。查找成果将回来基因模型称号（例如 evm.TU.Herato1801.64）及其方位（例如支架“Herato1801”坐落“1239943”到“1251211”方位）。
# optix H. erato
Herato1801:1239943-1251211
# optix H. melpomene
Hmel218003:705604-706407
当你点击基因模型时，你能够看到 optix 基因周围有哪些基因。您还能够在左侧看到“导出数据”按钮。这允许您将序列导出为 .fasta 文件。运用此功用，您不仅能够尝试导出 optix 基因，还能够导出它周围的 2,000,000 bp 区域。
您还能够在此处找到这些 .fasta 文件。
# scaffold Herato1801 start 1 end 2000000
sequences/Herato1801_1_2000000.fasta
# scaffold Hmel213003 start 1 end 2000000
sequences/Hmel213003_1_2000000.fasta
3. seq-seq-pan aligner
咱们能够运用 seq-seq-pan 将 fasta 文件中的序列组装成 pan 基因组。 Seq-seq-pan 经过构建复合共有序列或泛基因组来扩展多基因组比对器渐进式 Mauve 的功用，其间包括同源序列或部分共线块 (LCB) 以及每个基因组中的谱系特异性（非同源）序列基因组。然后将该泛基因组用作多基因组比对的参考坐标空间，其间包括任何基因组特有的序列。
关于咱们的序列，咱们将运用 seq-seq-pan，如下所示：
seq-seq-pan-wga --config genomefile=genome_list.txt outfilename=seq-seq-pan_out/SeqSeqPan_erato_melp_optix
Genome_list.txt 文件包括要包括在泛基因组组装中的 fasta 序列列表（每行一个）。该文件能够在这儿下载。
Seq-seq-pan 将输出几个文件。其间有两个与咱们相关：
- _consensus.fasta 文件包括共有泛基因组的完好 fasta 序列（将所有非同源序列拼接到组件中，并采用多个比对基因组中最常见的等位基因）。该一致文件划分了泛基因组的坐标空间，当咱们想要将原始基因组中的任何方位（例如TE方位）映射到泛基因组时将运用该一致文件。
- .xmfa 文件包括部分共线块 (LCB) 的列表。咱们将运用此文件来辨认同源或物种特异性的序列。
- 1：genome_list.txt 文件中第一个基因组的序列标识符。
- 2：genome_list.txt 文件中第二个基因组的序列标识符。
-
=区别单独的 LCB。 -
-对齐的 LCB 中存在间隙。
就是这样，咱们有了泛基因组！
4. 同享和共同的序列
咱们现在能够尝试确认序列的哪些部分在泛基因组中被辨认为同源或物种特异性。为此，咱们将运用自定义 Python 脚本，可在此处获取。
# For the Python script to work, we first need to remove newlines (enters) in the file from lines that include the sequence. This can be done with this perl oneliner:
perl -pe 'chomp if /^[ATCGNSBDHVMRWYK-]/' seq-seq-pan_out/SeqSeqPan_erato_melp_optix.xmfa| sed 's/=/n=/g' | sed 's/>/n>/g' | sed '/^$/d' > seq-seq-pan_out/SeqSeqPan_erato_melp_optix.noNewline.xmfa
# Now we can run the python script ('-g 1,2' is a list of the genome identifiers in the order of the genomes_list.txt file):
python seq-seq-pan_blocks_intervals.py -I seq-seq-pan_out/SeqSeqPan_erato_melp_optix.noNewline.xmfa -g 1,2
# The script annoyingly produces a 'TAB' at the end of each line which would break bedtools in the next step. We can remove that as follows:
sed -i 's/[[:space:]]*$//' seq-seq-pan_out/1_blocks_intervals.bed
sed -i 's/[[:space:]]*$//' seq-seq-pan_out/2_blocks_intervals.bed
输出是带有 | 的文件染色体|开端 |完毕 |每个基因组中序列的方位，但在泛基因组的坐标空间中（因而，当该序列被另一个基因组中的物种特异性序列打断时，会生成一条新线）。这种具有开端和完毕方位的格局通常称为 .bed 文件格局。
接下来，咱们能够运用bedtools 这些文件，得到每个基因组中的共同序列。
bedtools subtract -sorted -a 1_blocks_intervals.bed -b 2_blocks_intervals.bed > blocks_unique_1.bed
bedtools subtract -sorted -b 1_blocks_intervals.bed -a 2_blocks_intervals.bed > blocks_unique_2.bed
咱们还有一个自定义 Python 脚原本计算 LCB 内的序列同一性，可在此处获取。
# First, we need to transform the .xmfa file to a .fasta alignment
python seq-seq-pan_toFasta.py -I seq-seq-pan_out/SeqSeqPan_erato_melp_optix.noNewline.xmfa -g 1,2
cat genome* > seq-seq-pan_out/SeqSeqPan_erato_melp_optix.fasta
# Next, we need to find the sequences that are shared between genomes. We can do this with bedtools intersect.
bedtools intersect -sorted -a seq-seq-pan_out/1_blocks_intervals.bed -b seq-seq-pan_out/2_blocks_intervals.bed > seq-seq-pan_out/blocks_shared_1_2.bed
# Finally, we can calculate the conservation. This will output a .bed like file with identity scores for each shared interval.
python seq-seq-pan_bedfile_conservation.py -I seq-seq-pan_out/SeqSeqPan_erato_melp_optix.fasta -g 1,2 -b seq-seq-pan_out/blocks_shared_1_2.bed -o seq-seq-pan_out/conservation_shared_1_2.bed
5. 将注释映射到泛基因组
seq-seq-pan 的映射功用允许将所包括基因组的任何原始方位转换为泛基因组（=泛基因组坐标）。该函数将一个文件作为输入，该文件包括单列方位和第一行，该文件指定从何处映射到何处（例如 2tc，这意味着从基因组 2 进行映射（Hmel218003 序列，它是基因组列表中的第二个基因组） .txt 文件）到基因组 c（一致泛基因组序列））。
这儿我有 Hmel218003 支架上 optix 基因的开始和完毕方位（您也能够从 ensembl.lepbase.org/ 获取这些或任何基因的方位。
2 c
705604
706407
运转映射函数：
seq-seq-pan map -c seq-seq-pan_out/SeqSeqPan_erato_melp_optix_consensus.fasta -p ./ -i gene_annotation/optix_Hmel2.toMap.txt -n gene_annotation/optix_Hmel2.toMap.pan
现在让咱们在 R 中绘制泛基因组。
本文由mdnice多渠道发布
