6
$\begingroup$

I have this GTF file and I use the command below on a Linux machine to extract the coding regions of the genome:

awk '{if($3=="transcript" && $20=="\"protein_coding\";"){print $0}}' gencode.gtf

How I could do the inverse and keep only non coding regions?

$\endgroup$
  • 2
    $\begingroup$ Do you want all non-coding regions of the genome or do you want all non-coding transcripts? These are two very different things. $\endgroup$ – terdon Feb 22 '19 at 18:24
4
$\begingroup$

Getting the non coding regions of a protein coding transcript, sounds like you are looking for UTR.

UTR has its own feature in the gtf file. So you can do this:

$ awk -v FS="\t" '$3=="UTR"' gencode.gtf

If the gtf file is compressed use this instead:

$ zcat gencode.gtf.gz | awk -v FS="\t" '$3=="UTR"'

BTW: Why are you using such an old release of gencode? The current version is v29.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Sorry, literally I need non coding regions of human genome, but for asking my question here I referred to coding parts too $\endgroup$ – Exhausted Feb 22 '19 at 12:53
  • $\begingroup$ Sorry I tried hat but my output is empty $\endgroup$ – Exhausted Feb 22 '19 at 12:59
  • 1
    $\begingroup$ As @Wouter tells you, the non coding region of a genome is the complement of the coding regions. Coding regions have its own feature in the gtf file. You can get them with $ awk -v FS="\t" '$3=="CDS"' gencode.gtf. Reading the manual for bedtools complement is your task. $\endgroup$ – finswimmer Feb 22 '19 at 13:14
  • $\begingroup$ Sorry but your commands return nothing, I mean not working returning empty file $\endgroup$ – Exhausted Feb 22 '19 at 18:24
  • $\begingroup$ The gtf file the OP has linked to includes non-coding transcripts (LINCs, pseudogenes, tRNAs etc). I am guessing this is what they're after. $\endgroup$ – terdon Feb 22 '19 at 18:37
2
$\begingroup$

This isn't a problem that's easily solved with awk. It's not like you're extracting a feature that's annotated in the GTF file. Instead, you want the empty space between annotated features.

A few years ago I wrote a program called LocusPocus for a similar task. It uses a gene annotation to break down a genome into gene loci and intergenic regions. It handles overlapping annotations and other weirdness pretty robustly. The output will include both coding regions and non-coding regions, but you can identify the intergenic spaces as those with iLocus_type equal to iiLocus or fiLocus.

Note: the --delta parameter will extend each gene/transcript by 500bp by default.

Caveat: the program only accepts GFF3 input by default. Hopefully it won't be too hard to convert your GTF to GFF3.

Another caveat: eventual interpretation of these data will depend on what features are annotated in the genome and which annotations you include vs ignore. Do you want your non-coding regions to include non-coding genes, or should these be treated separately? Some non-coding regions will be full of transposable elements and other repetitive DNA, while others will have enhancers, promoters, or other regulatory elements. It's important to tread carefully before you jump to any conclusions.

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ Absolutely brilliant name! :) $\endgroup$ – terdon Feb 22 '19 at 19:32
  • $\begingroup$ Sorry, thank you for explanation; After googling I see I actually need enhancers, promoters, or other regulatory elements of human genome to find cancer driver genes placed in these regions. $\endgroup$ – Exhausted Feb 24 '19 at 21:24
2
$\begingroup$

If you want all transcripts from that gtf file whose type isn't "protein_coding", you can use almost the same command, just change the == ("is") to != ("isn't"):

awk '{if($3=="transcript" && $20!="\"protein_coding\";"){print $0}}' gencode.gtf 

Or, a simpler version:

awk '$3=="transcript" && $20!="\"protein_coding\";"' gencode.gtf 

Note that this will not include any of the havana transcripts in the file, but I am assuming that's what you want since that's what your original command did.

Specifically, the command will return the following types of transcript (the numbers on the left are the number of such transcripts in the file):

awk '$3=="transcript" && $20!="\"protein_coding\";"{print $20}' gencode.gtf  | sort | uniq -c | sort -nk1
      1 "translated_processed_pseudogene";
      2 "Mt_rRNA";
      3 "IG_J_pseudogene";
      3 "TR_D_gene";
      4 "TR_J_pseudogene";
      5 "TR_C_gene";
     10 "IG_C_pseudogene";
     18 "IG_C_gene";
     18 "IG_J_gene";
     22 "Mt_tRNA";
     25 "3prime_overlapping_ncrna";
     27 "TR_V_pseudogene";
     37 "IG_D_gene";
     58 "non_stop_decay";
     59 "polymorphic_pseudogene";
     74 "TR_J_gene";
     97 "TR_V_gene";
    144 "IG_V_gene";
    182 "unitary_pseudogene";
    196 "IG_V_pseudogene";
    330 "sense_overlapping";
    387 "pseudogene";
    442 "transcribed_processed_pseudogene";
    531 "rRNA";
    802 "sense_intronic";
    860 "transcribed_unprocessed_pseudogene";
   1529 "snoRNA";
   1923 "snRNA";
   2050 "misc_RNA";
   2549 "unprocessed_pseudogene";
   3116 "miRNA";
   9710 "antisense";
  10623 "processed_pseudogene";
  11780 "lincRNA";
  13052 "nonsense_mediated_decay";
  25955 "retained_intron";
  28082 "processed_transcript";

You might also want to remove that "translated_processed_pseudogene" since that is actually translated into protein and is therefore technically coding:

awk '$3=="transcript" && 
     $20!="\"protein_coding\";" &&
     $20!="\"translated_processed_pseudogene\";"' gencode.gtf  
|improve this answer|||||
$\endgroup$
  • $\begingroup$ Thanks a lot, really thank you for saving me I could not solve that myself. How I can extract the below information from each line of resulting non-coding file chr1 29553 30039 ENSG00000243485.2 + gene_name "MIR1302-11" $\endgroup$ – Exhausted Feb 22 '19 at 20:29
  • $\begingroup$ @FereshTeh you're welcome. I think you want awk '$3=="transcript" && $20!="\"protein_coding\";" && $20!="\"translated_processed_pseudogene\";"{print $1,$4,$5,$10,$7}' gencode.gtf but, if not, please ask a new question about that. $\endgroup$ – terdon Feb 22 '19 at 20:32
  • $\begingroup$ Thanks a lot that returns all except gene name, this is output chr1 29554 31097 "ENSG00000243485.2"; + $\endgroup$ – Exhausted Feb 22 '19 at 21:15
  • 1
    $\begingroup$ @FereshTeh please ask a new question so you can show exactly what output you need. $\endgroup$ – terdon Feb 22 '19 at 21:17
0
$\begingroup$

Likely non-coding regions of genome are here

https://elifesciences.org/download/aHR0cHM6Ly9jZG4uZWxpZmVzY2llbmNlcy5vcmcvYXJ0aWNsZXMvMjE3NzgvZWxpZmUtMjE3Nzgtc3VwcDQtdjMueGxzeA==/elife-21778-supp4-v3.xlsx?_hash=KQi5jfO3kT2c4Qw44j4Rg6YAyCBQilYuWHVYXcRDuuo%3D

|improve this answer|||||
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.