Sollevamento da hg38 a hg37

Aug 26 2020

Sto riscontrando strani problemi nel sollevare questo set di dati da hg37 a hg38. Sto usando i dati HGDP WGS da ftp://ngs.sanger.ac.uk/production/hgdp/hgdp_wgs.20190516/ . Il set di dati è stato diviso per cromosoma, quindi li ho prima convertiti in BED e poi li ho uniti in un singolo file usando plink1.9. Tuttavia, quando provo a passare a hg37/19 utilizzando liftOver o CrossMap v0.5.1 ottengo strani errori. Ho provato sia la catena da ensemble.org che da UCSC , le ho provate sia compresse che non compresse. Comincio a pensare di dover provare la conversione sui file VCF ma questi sono così grandi che preferirei non duplicarli inutilmente.

Dati

I dati sembrano a posto. Se guardiamo al BIM:

user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38_chrpos.bim
1   1:10153 0   10153   G   A
1   1:10163 0   10163   C   T
1   1:10180 0   10180   C   T
1   1:10250 0   10250   C   A
1   1:10257 0   10257   C   A
1   1:10291 0   10291   T   C
1   1:10297 0   10297   T   C
1   1:10327 0   10327   C   T
1   1:10330 0   10330   A   C
1   1:10333 0   10333   T   C
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail hgdp_hg38_chrpos.bim
Y   Y:56887228  0   56887228    C   T
Y   Y:56887285  0   56887285    C   A
Y   Y:56887316  0   56887316    A   G
Y   Y:56887463  0   56887463    A   C
Y   Y:56887491  0   56887491    A   C
Y   Y:56887503  0   56887503    T   G
Y   Y:56887583  0   56887583    G   A
Y   Y:56887631  0   56887631    T   C
Y   Y:56887837  0   56887837    A   G
Y   Y:56887844  0   56887844    C   T

I dati sono nei nomi delle varianti chrpos qui, ma ho avuto lo stesso risultato usando il rsid originale.

Catene

Le catene sembrano a posto e sono in grado di leggerle in R senza problemi:

user@desktop:/media/luks8tb1/data/genomics/HGDP$ head GRCh38_to_GRCh37.chain
chain 1 1 248956422 + 10000 297968 1 249250621 + 10000 267719 2
167417 80249 50000
40302

chain 1 1 248956422 + 347968 501617 1 249250621 - 248779253 248932902 3
153649

chain 1 1 248956422 + 585988 12949384 1 249250621 + 521368 13009210 4
1044707 1 0
3709 3 0
user@desktop:/media/luks8tb1/data/genomics/HGDP$ tail GRCh38_to_GRCh37.chain

chain 1 chrY 57227415 + 56821509 57217415 Y 59373566 + 58967656 59363566 11183
85168 0 1
78580 0 1
63496 0 1
3158 1 0
30382 0 1
54147 0 1
80974

CrossMap

user@desktop:/media/luks8tb1/data/genomics/HGDP$ CrossMap.py bed GRCh38_to_GRCh37.chain hgdp_hg38_chrpos.bed hgdp_hg37_chrpos.bed
@ 2020-08-25 20:04:26: Read the chain file:  GRCh38_to_GRCh37.chain
Traceback (most recent call last):
  File "/home/user/.local/bin/CrossMap.py", line 166, in <module>
    crossmap_bed_file(mapTree, in_file, out_file)
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/mapbed.py", line 32, in crossmap_bed_file
    for line in ireader.reader(inbed):
  File "/home/user/.local/lib/python3.6/site-packages/cmmodule/ireader.py", line 27, in reader
    yield l.decode('utf8').strip().replace("\r", "")
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xdd in position 20: invalid continuation byte

liftOver

user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.bed hg38ToHg19.over.chain hgdp_hg37.bed hgdp_hg37_errors.bed
Reading liftover chains
Mapping coordinates
Data format error: expecting at least 3 fields in BED file (hgdp_hg38.bed)

Provando il file PED con liftOver

Ho convertito il file BED in PED con plink1.9:

plink --bfile hgdp_hg38 --recode --out hgdp_hg38

Quindi:

user@desktop:/media/luks8tb1/data/genomics/HGDP$ liftOver hgdp_hg38.ped hg38ToHg19.over.chain hgdp_hg37.ped hgdp_hg37_errors.ped
Reading liftover chains
Mapping coordinates
invalid unsigned integer: "HGDP00001"

Ci sono alcune domande su Google a riguardo, ma niente che sembri utile.

CrossMap non supporta PED.

Risposte

1 PhoenixMu Aug 27 2020 at 20:39

Il formato BED può fare riferimento al formato di file ".txt" che memorizza i cromosomi e le posizioni delle caratteristiche genomiche, nonché alcune altre cose come il nome. bedformato può riferirsi anche al formato plink, che appartiene a bed, bim, fam trio. Il formato plink bed memorizza tutte le varianti genetiche, ma il programma deve fare riferimento al file bim per le loro posizioni.

Ovviamente il liftover del genoma richiede solo il coordinamento di ciascun SNP e non necessita affatto di singoli genotipi. Infatti, il formato del letto in PLINK non fornisce alcuna informazione utile. Sapere questo significa che fornire il file del letto plink è sbagliato. E questo è abbastanza chiaro dal messaggio di errore.

Inoltre, poiché i tre file plink devono sempre essere usati insieme, la maggior parte dei programmi li chiama formato file plink . Quindi ogni volta che le persone dicono formato LETTO, è il formato del file txt.

ATpoint Aug 26 2020 at 04:45

Non sono sicuro di crossMap ma liftOver richiede il formato BED dove la colonna 1 è il cromosoma, la colonna 2 l'inizio basato su 0 e la colonna 3 la coordinata finale. Le colonne aggiuntive sono facoltative e possono essere utilizzate per memorizzare le altre colonne che hai. Basta convertire BIM in BED ed eseguire liftOver.