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
1 Answer
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
-
$\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$– bliApr 11, 2019 at 12:06
p13
refers to in kayvee's answer below. $\endgroup$