Add CHESS to IGB Quickload for hg38

Contacted authors to find out location of reference assembly used to generate the annotations. Answer is:

Hi Ann and Steven,

We used CHM13 2.0 as our target sequence, which you can find here: https://github.com/marbl/CHM13.
All chromosomes except chromosome Y in CHM13 v1.1 are the same as CHM13 v2.0. ("Changes from v1.1 is the addition of a finished chromosome Y from the GIAB HG002/NA24385 sample, sequenced both by GIAB and HPRC.")

Best regards,
Kuan-Hao

Next step:

• locate candidate assembly sequence we could deploy to an igbquickload site that would enable IGB users to view the CHESS gene model annotations

Locating candidate CH13 genome sequence:

• This is "T2T-CHM13v2.0" https://github.com/marbl/CHM13#t2t-chm13v20
• Retrieved genome sequence data file with:

wget https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0.fa.gz 
--2023-03-03 16:19:58--  https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/analysis_set/chm13v2.0.fa.gz
Resolving s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)... 52.218.228.208, 52.218.176.72, 52.218.229.200, ...
Connecting to s3-us-west-2.amazonaws.com (s3-us-west-2.amazonaws.com)|52.218.228.208|:443... connected.
HTTP request sent, awaiting response... 200 OK
Length: 981185361 (936M) [binary/octet-stream]
Saving to: ‘chm13v2.0.fa.gz’

Notes on size differences:

• Unpacked 0.9 Gb chm13v2.0.fa.gz to 2.9 Gb chm13v2.0.fa
• Converted 2.9 Gb chm13v2.0.fa to 0.90 Gb chm13v2.0.2bit

UCSC has this assembly, too. UCSC long name for CH13v2.0 is "human T2T-CHM13 v2.0 (GCA_009914755.4)". This long name belongs in the genome assembly description field of the quickload contents table.

Start using variety/cultivar trinomial to indicate assembly.

For example:

• H_sapiens_CHM13v2_Mar_2022

Ref: The complete sequence of a human genome

CHM13v2 2bit file: chess/H_sapiens_CHM13v2_Mar_2022/chm13v2.0.2bit deployed to RENCI host.

Downloaded GFF file for hg38 and opened in IGB. Looks OK.
Checked column 2. Values included:

gunzip -c chess3.0.gff.gz | cut -f2 | sort | uniq
BestRefSeq
Curated Genomic
ENSEMBL
FANTOM
Gnomon
HAVANA
RefSeq
StringTie
cmsearch
ensembl
ensembl_havana
havana
tRNAscan-SE

Why HAVANA and havana? Upper-case and lower-case both present - why?

Sample transcript information:

chr1    BestRefSeq      transcript      34611   36081   .       -       .       ID=CHS.4.1;geneID=CHS.4;gene_type=lncRNA;gene_name=FAM138A;db_xref=RefSeq:NR_026818.1,GENCODE:ENST00000417324.1;num_samples=233;max_tpm=6.158023;assembly_id=ALL_00180370

Looks like there's a lot of useful computational results included in this extra feature field.
Normally I put the gene description from NCBI Entrez Gene into field 14 of the BED-detail data file for IGB, but probably I don't need to do this here because the reference annotations from UCSC are already loaded. Probably it would be better to instead use this entire "extra feature" field instead, except geneID, which I can insert into field 13, where I normally insert gene symbol(s).

Used genomesource code to create a list of compound dna seqfeature objects. Exploring:

>>> model.getProducer()
>>> model.getFeatType()
'mRNA'
>>> model.getKeyVals()
{'producer': 'BestRefSeq', 'db_xref': 'RefSeq:NM_001174081.3_1,GENCODE:ENST00000621200.4', 'geneID': 'CHS.27883', 'gene_type': 'protein_coding', 'ID': 'CHS.27883.1', 'gene_name': 'NLRP2'}
>>> model.getVal("Parent")
>>> model.getVal("geneID")
'CHS.27883'
>>> geneIDs = map(lambda x:x.getVal("geneID"),models)
>>> len(geneIDs)
168475
>>> len(models)
168475
>>> map(lambda x: x.getFeatType(),model.getFeats())
['exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'exon', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS', 'CDS']

The number of models created matches the number of "transcript" features in the data file

Made a first-pass BED14 file from hg38 CHESS gene model annotations. Checked sequence per assembly number before and after loading test file. It was 595 before and after. Same with loading original CHESS GFF file.

In CHESS, RAP1GAP on chromosome 1 has lots of splice variants. The NCBI Gene record (abstract is below) mentions expression and mis-regulation in diverse cancer types. Also, looks like the UCSC reference annotations IGB shows probably need to be updated, as IGB does not show nearly as many gene models as visible and mentioned in the Entrez Gene record.

This gene encodes a type of GTPase-activating-protein (GAP) that down-regulates the activity of the ras-related RAP1 protein. RAP1 acts as a molecular switch by cycling between an inactive GDP-bound form and an active GTP-bound form. The product of this gene, RAP1GAP, promotes the hydrolysis of bound GTP and hence returns RAP1 to the inactive state whereas other proteins, guanine nucleotide exchange factors (GEFs), act as RAP1 activators by facilitating the conversion of RAP1 from the GDP- to the GTP-bound form. In general, ras subfamily proteins, such as RAP1, play key roles in receptor-linked signaling pathways that control cell growth and differentiation. RAP1 plays a role in diverse processes such as cell proliferation, adhesion, differentiation, and embryogenesis. Alternative splicing results in multiple transcript variants encoding distinct proteins. [provided by RefSeq, Aug 2011]

Committed new version of GFF3 / BED parsing / writing code.

Repository: https://bitbucket.org/lorainelab/genomesource/src/master/
Commit: https://bitbucket.org/lorainelab/genomesource/commits/92eb0500f7cefd7b86ecbefa3d763a4d1e5f5c34

Made BED14 file as follows:

local aloraine$ CHESSgff3ToBedDetail.py -g chess3.0.gff.gz -b chess3.0.GRCh38.p12.bed 
local aloraine$ wc -l chess3.0.GRCh38.p12.bed 
  168475 chess3.0.GRCh38.p12.bed
local aloraine$ sort -k1,1 -k2,2n chess3.0.GRCh38.p12.bed | bgzip > chess3.0.GRCh38.p12.bed.gz 
local aloraine$ tabix -s 1 -b 2 -e 3 chess3.0.GRCh38.p12.bed.gz
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
[W::tbx_parse1] Coordinate <= 0 detected. Did you forget to use the -0 option?
local aloraine$ tabix -0 -s 1 -b 2 -e 3 -f chess3.0.GRCh38.p12.bed.gz

Added chess 3.0 BED-detail file to svn repository, updated to revision 170.
See: https://svn.bioviz.org/viewvc/genomes/quickload/H_sapiens_Dec_2013/

Updated igbquickload at RENCI:

[aloraine@lorainelab-quickload quickload]$ svn up
Updating '.':
A    H_sapiens_Dec_2013/chess3.0.GRCh38.p12.bed.gz
U    H_sapiens_Dec_2013/annots.xml
A    H_sapiens_Dec_2013/chess3.0.GRCh38.p12.bed.gz.tbi
Updated to revision 170.

See: http://lorainelab-quickload.scidas.org/quickload/H_sapiens_Dec_2013/

I'm pretty happy with how the gene models look in IGB. I am going to move the work needed to set up CHESS for CHM13 into a new ticket linked with this one.

Wrote and submitted review. I liked the article and recommended to "accept". I also mentioned my confusion about the GFF producer and extra feature fields, and suggested adding more documentation to the Web site.