In [1]:
%run ~/relmapping/annot/notebooks/annot__init__.ipynb
annot_ = 'annot_ce11'
def mp(fp, annot_=annot_): return os.path.join(annot_, 'canonical_geneset', fp)
In [2]:
# Read WB260 canonical geneset
kwargs_to_gtf = {'index': False, 'header': False, 'quoting': csv.QUOTE_NONE, 'sep': '\t'}
fp_ = 'wget/ftp.wormbase.org/pub/wormbase/releases/WS260/species/c_elegans/PRJNA13758/c_elegans.PRJNA13758.WS260.canonical_geneset.gtf.gz'
df_ = pd.read_csv(fp_, sep='\t', names=yp.NAMES_GTF, comment='#')
# Strip mitochondrial & change chromosome names to ce11
df_ = df_.query('(chrom != "CHROMOSOME_MtDNA") & (chrom != "MtDNA")').sort_values(['chrom', 'start', 'end', 'strand']).reset_index(drop=True)
df_['chrom'] = list(map(lambda chrom: 'chr' + chrom.lstrip('CHROMOSOME_'), df_['chrom']))
df_.head()
Out[2]:
In [3]:
# locus_id: intended to roughly match labels as shown on the default IGV RefSeq track
def locus_id_from_gene_id(l_gene_id):
# Attach display_id
def locus_id_(locus, sequence_id, gene_id):
if locus == locus:
return locus
elif sequence_id == sequence_id:
return sequence_id
else:
return gene_id
fp_geneIDs = 'wget/ftp.wormbase.org/pub/wormbase/releases/WS260/species/c_elegans/PRJNA13758/annotation/c_elegans.PRJNA13758.WS260.geneIDs.txt.gz'
l_ = ['gene_id', 'locus', 'sequence_id', 'status']
df_geneIDs = pd.read_csv(fp_geneIDs, sep=',', names=('na', 'gene_id', 'locus', 'sequence_id', 'status'))[l_]
df_geneIDs['display_id'] = list(map(locus_id_,df_geneIDs['locus'], df_geneIDs['sequence_id'], df_geneIDs['gene_id']))
df_geneIDs = df_geneIDs.set_index('gene_id')
for gene_id in l_gene_id:
yield df_geneIDs.loc[gene_id]['display_id']
# Parse attribute string
df_attr_ = pd.DataFrame(df_['attribute'].apply(hts.parse_GFF_attribute_string).tolist())
# Define locus_id
se_locus_id = [*locus_id_from_gene_id(df_attr_['gene_id'])]
# Attach locus_id at the end of the unparsed attr field
def fmt_(attribute, locus_id): return '%s locus_id "%s";' % (attribute, locus_id)
df_['attribute'] = list(map(fmt_, df_['attribute'], se_locus_id))
df_.head()
Out[3]:
In [ ]:
df_['feature'].value_counts()
Out[ ]:
In [ ]:
# WS260_ce10.canonical_geneset.gtf.gz: entire canonical geneset
yp.to_csv_gz(mp('WS260_ce11.canonical_geneset.gtf.gz'), df_, **yp.TO_GTF_KWARGS)
!gunzip -c {mp('WS260_ce11.canonical_geneset.gtf.gz')} | wc -l
In [ ]:
# WS260_ce10.genes.gtf.gz: all gene records
yp.to_csv_gz(mp('WS260_ce11.genes.gtf.gz'), df_.query('feature=="gene"'), **yp.TO_GTF_KWARGS)
!gunzip -c {mp('WS260_ce11.genes.gtf.gz')} | wc -l
In [ ]:
# WS260_ce10.transcripts.gtf.gz: all transcript records
yp.to_csv_gz(mp('WS260_ce11.transcripts.gtf.gz'), df_.query('feature!="gene"'), **yp.TO_GTF_KWARGS)
!gunzip -c {mp('WS260_ce11.transcripts.gtf.gz')} | wc -l
In [ ]:
df_attr_['transcript_biotype'].value_counts()
In [ ]:
# WS260_ce11.annot.gtf.gz: all transcripts used in the annotation
l_annot_ = ['protein_coding', 'pseudogene', 'tRNA', 'snoRNA', 'miRNA', 'snRNA', 'rRNA']
m_annot_ = df_attr_['transcript_biotype'].isin(l_annot_)
fp_ = os.path.join(annot_, 'WS260_ce11.annot.gtf.gz')
yp.to_csv_gz(fp_, df_[m_annot_].query('feature!="gene"'), **yp.TO_GTF_KWARGS)
!gunzip -c {fp_} | wc -l