Для анализа крупных эволюционных событий были взяты геномы 5 штаммов Bradyrhizobium japonicum. Ниже приведен файл genomes.tsv в котором указаны: база данных, AC, штамм, хромосома, тип ДНК (кольцевая/линейная) и полное название штамма.
! cat src/genomes.tsv
all:embl:CP058354.1 5038 chr1 c Bradyrhizobium japonicum strain 5038 all:embl:AP012206.1 USDA-6 chr1 c Bradyrhizobium japonicum strain USDA 6 all:embl:CP017637.1 J5 chr1 c Bradyrhizobium japonicum strain J5 all:embl:CP010313.1 E109 chr1 c Bradyrhizobium japonicum strain E109 all:embl:CP007569.1 SEMIA-5079 chr1 c Bradyrhizobium japonicum strain SEMIA 5079
Постройка нуклеотидного пангенома происходила по последовательности команд ниже.
C:\Users\jakew\Desktop\npg\JW>npge -g npge.conf
C:\Users\jakew\Desktop\npg\JW>npge Prepare > log_prep.txt 2>&1
C:\Users\jakew\Desktop\npg\JW>npge Examine > log_ex.txt 2>&1
C:\Users\jakew\Desktop\npg\JW>echo "MIN_IDENTITY = 0.864"
"MIN_IDENTITY = 0.864"
C:\Users\jakew\Desktop\npg\JW>npge MakePangenome > log_pg.txt 2>&1
C:\Users\jakew\Desktop\npg\JW>npge PostProcessing > log_pp.txt 2>&1
C:\Users\jakew\Desktop\npg\JW>qnpgeЛоги и конфиг после расчетов:
Основные файлы выдачи:
Для решения данной задачи использовался язык программирования Python и библиотека Pandas.
import pandas as pd
data = pd.read_csv('src/pangenome.bi', sep = '\t')
Ниже представлены первые 5 блоков, которые храняться в файле pangenome.bi.
data.head(5)
| block | fragments | cols | ident-nogap | ident-gap | noident-nogap | noident-gap | pure-gap | ident | GC | 5038 | E109 | J5 | SEMIA-5079 | USDA-6 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | r55x462 | 55 | 462 | 445 | 6 | 9 | 2 | 0 | 0.9696 | 0.5933 | 8 | 8 | 15 | 14 | 10 |
| 1 | r54x272 | 54 | 272 | 271 | 0 | 1 | 0 | 0 | 0.9963 | 0.5735 | 8 | 8 | 14 | 14 | 10 |
| 2 | r54x105 | 54 | 105 | 105 | 0 | 0 | 0 | 0 | 1.0000 | 0.6095 | 8 | 8 | 14 | 14 | 10 |
| 3 | r48x338 | 48 | 338 | 317 | 17 | 4 | 0 | 0 | 0.9630 | 0.5505 | 9 | 9 | 10 | 11 | 9 |
| 4 | r47x262 | 47 | 262 | 256 | 4 | 2 | 0 | 0 | 0.9847 | 0.6261 | 9 | 9 | 9 | 11 | 9 |
data.info()
<class 'pandas.core.frame.DataFrame'> RangeIndex: 9407 entries, 0 to 9406 Data columns (total 15 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 block 9407 non-null object 1 fragments 9407 non-null int64 2 cols 9407 non-null int64 3 ident-nogap 9407 non-null int64 4 ident-gap 9407 non-null int64 5 noident-nogap 9407 non-null int64 6 noident-gap 9407 non-null int64 7 pure-gap 9407 non-null int64 8 ident 9407 non-null float64 9 GC 9407 non-null float64 10 5038 9407 non-null int64 11 E109 9407 non-null int64 12 J5 9407 non-null int64 13 SEMIA-5079 9407 non-null int64 14 USDA-6 9407 non-null int64 dtypes: float64(2), int64(12), object(1) memory usage: 1.1+ MB
Для поиска самых длинных делеций сначала были убраны ненужные переменные, а из наблюдений оставлены только h-блоки.
h_blocks = data.loc[data.block.str.contains('h')].sort_values('cols', ascending = 0).iloc[:, [0, 2, 10, 11, 12, 13, 14]]
h_blocks
| block | cols | 5038 | E109 | J5 | SEMIA-5079 | USDA-6 | |
|---|---|---|---|---|---|---|---|
| 5292 | h3x42062 | 42062 | 1 | 1 | 0 | 0 | 1 |
| 5293 | h3x35534 | 35534 | 1 | 1 | 0 | 0 | 1 |
| 6004 | h2x32417 | 32417 | 0 | 0 | 1 | 1 | 0 |
| 6005 | h2x25936 | 25936 | 0 | 0 | 1 | 1 | 0 |
| 5294 | h3x21125 | 21125 | 1 | 1 | 0 | 0 | 1 |
| ... | ... | ... | ... | ... | ... | ... | ... |
| 5794 | h3x100n2 | 100 | 1 | 1 | 0 | 0 | 1 |
| 5795 | h3x100n3 | 100 | 1 | 1 | 0 | 0 | 1 |
| 5796 | h3x100n4 | 100 | 1 | 1 | 0 | 0 | 1 |
| 5797 | h3x100n5 | 100 | 1 | 1 | 0 | 0 | 1 |
| 6465 | h2x100n1 | 100 | 0 | 0 | 1 | 1 | 0 |
2117 rows × 7 columns
После была создана новая переменная sum, которая считает сколько штаммов имеют данный блок у себя в геноме. Так как нужно найти такие блоки, чтобы их не было только в одном штамме, данные были отфильтрованы по значению переменной sum равному 4, то есть данный блок присутвует у 4 штаммов, а у одного - нет.
Сначала производился поиск делеций в штамме 5038, однако, по результатам фильтрации, оказалось, что данный штамм никогда не терял блок в одиночку. Всегда есть хотя бы один штамм, который терял этот же блок.
h_blocks['sum'] = h_blocks.apply(lambda row: row['5038' : 'USDA-6'].sum(), axis=1)
h_blocks.loc[(h_blocks['sum'] == 4) & (h_blocks['5038'] == 0)]
| block | cols | 5038 | E109 | J5 | SEMIA-5079 | USDA-6 | sum |
|---|
Далее штамм E109. Самой длинной делецией для данного штамма оказалась потеря блока h4x227n3 длинной 227. В этом блоке находились гены: у 5038 нет аннотации, BKD09_47320 transposase (J5), BJS_08812 transposase (SEMIA-5079) и BJ6T_88600 transposase (USDA-6).
h_blocks.loc[(h_blocks['sum'] == 4) & (h_blocks['E109'] == 0)].iloc[0, :]
block h4x227n3 cols 227 5038 1 E109 0 J5 1 SEMIA-5079 1 USDA-6 1 sum 4 Name: 3724, dtype: object
Штамм J5 потерял блок длиной 16381 h4x16381. Гены: нет аннотации (5038), RN69_01300 transcriptional regulator (E109), BJS_05560 hypothetical protein (SEMIA-5079), BJ6T_02700 hypothetical protein (USDA-6).
h_blocks.loc[(h_blocks['sum'] == 4) & (h_blocks['J5'] == 0)].iloc[0, :]
block h4x16381 cols 16381 5038 1 E109 1 J5 0 SEMIA-5079 1 USDA-6 1 sum 4 Name: 3115, dtype: object
Следующий штамм - SEMIA-5079 потерял блок длинной 9710 h4x9710. Гены: нет аннотации (5038), RN69_07820 3-oxoacyl-ACP reductase (E109), BKD09_11055 3-oxoacyl-ACP reductase (J5), BJ6T_15830 hypothetical protein (USDA-6).
h_blocks.loc[(h_blocks['sum'] == 4) & (h_blocks['SEMIA-5079'] == 0)].iloc[0, :]
block h4x9710 cols 9710 5038 1 E109 1 J5 1 SEMIA-5079 0 USDA-6 1 sum 4 Name: 3117, dtype: object
У последнего штамма произошла делеция длинной 8873 в блоке h4x8873. Гены: нет аннотации (5038), RN69_43045 conjugal transfer protein TrbB (E109), BKD09_47710 P-type conjugative transfer ATPase TrbB (J5), BJS_08460 Conjugal transfer protein (SEMIA-5079).
h_blocks.loc[(h_blocks['sum'] == 4) & (h_blocks['USDA-6'] == 0)].iloc[0, :]
block h4x8873 cols 8873 5038 1 E109 1 J5 1 SEMIA-5079 1 USDA-6 0 sum 4 Name: 3118, dtype: object
У штамма J5 во время эволюционного процесса произошла перестановка длока g5x20648 на 188 позицию, в то время как у остальных штаммов этот блок находится на 194 позиции.
Поиск производился по стабильным блокам. В блоке S5x12321 для штамма SEMIA-5079 указано, что эта последовательность Putative C4-dicarboxylate, в то время как у J5 и E109 - ABC transporter, а у остальных либо ничего (5038), либо гипотетический белок (USDA-6), при совпадении 95.71%.