1
$\begingroup$

I have been given a bunch of .bam files, that come from the alignment on GRCh37d5 of human genome by Star tool. I want to extract raw read counts from these bam files by FEATURECOUNTS tool so I need a GTF. I googled a lot but I did not find GTF for GRCh37d5. I used another GTFs and I obtained weird features. Please one of you help me in finding this specific GTF. I converted a bam to fastq and I used Salmon to quantify gene expression by a transcripome.fasta but rate of mapping was 0.06

$\endgroup$
  • 1
    $\begingroup$ It's helpful to know that GRCh37 refers to the Genome Reference Consortium assembly of the human genome, version 37 (should be equivalent to hg19 from UCSC). I've never seen "5d" though. Often times GRC will introduce "patch" updates with minor changes that don't change the primary coordinate system—that's what p13 refers to in kayvee's answer below. $\endgroup$ – Daniel Standage Apr 10 '19 at 18:30
  • 1
    $\begingroup$ @DanielStandage Might be of interest to you: biostars.org/p/328824/#330791 $\endgroup$ – Ram RS Apr 10 '19 at 18:58
2
$\begingroup$

The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo_sapiens/GRCh37.p13_interim_annotation/

You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:

$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses chr# format and Ensembl uses #. For this, you can use the -A parameter from featureCounts. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:

## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv
| improve this answer | |
$\endgroup$
  • $\begingroup$ I didn't know about "SAF" format. I use gff with FeatureCounts. For the record, according to rdocumentation.org/packages/Rsubread/versions/1.22.2/topics/…: "The SAF annotation format has five required columns, including GeneID, Chr, Start, End and Strand. These columns can be in any order". $\endgroup$ – bli Apr 11 '19 at 12:06

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.