Error with Rseqc infer_experiment - unknown datatype
I'm using rseqc on my samples bam file to check the strand specific. For this I used gencode gtf file. Initially I converted gencode gtf to bed file like:
This is how the gtf looks:
##description: evidence-based annotation of the human genome (GRCh38), version 27 (Ensembl 90)
##provider: GENCODE
##contact: gencode-help@sanger.ac.uk
##format: gtf
##date: 2017-08-01
chr1 HAVANA gene 11869 14409 . + . gene_id "ENSG00000223972.5"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; level 2; havana_gene "OTTHUMG00000000961.2";
chr1 HAVANA transcript 11869 14409 . + . gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; level 2; transcript_support_level "1"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1 HAVANA exon 11869 12227 . + . gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 1; exon_id "ENSE00002234944.1"; level 2; transcript_support_level "1"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1 HAVANA exon 12613 12721 . + . gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 2; exon_id "ENSE00003582793.1"; level 2; transcript_support_level "1"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
chr1 HAVANA exon 13221 14409 . + . gene_id "ENSG00000223972.5"; transcript_id "ENST00000456328.2"; gene_type "transcribed_unprocessed_pseudogene"; gene_name "DDX11L1"; transcript_type "processed_transcript"; transcript_name "DDX11L1-202"; exon_number 3; exon_id "ENSE00002312635.1"; level 2; transcript_support_level "1"; tag "basic"; havana_gene "OTTHUMG00000000961.2"; havana_transcript "OTTHUMT00000362751.1";
I created bed file with:
awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' gencode.v27.annotation.gtf | gtf2bed - > gencode.v27.annotation.bed
And then I used that bed file for infer experiment.py
python infer_experiment.py -r gencode.v27.annotation.bed -I Sample.sorted.bam
Error:
Reading reference gene model gencode.v27.annotation.bed ... Done
Loading SAM/BAM file ... Finished
Total 0 usable reads were sampled
Unknown Data type
1 answer
Log in to answer this question.
More posts like this
-
bioawk command fails to generate tx2gene.tsv
written by Claire •
-
Is there a way to safely remove all underscore from gene id for hg19 from genecode
written by simplitia •
-
Counting how many transcripts per gene in gtf file
written by sallyey2 •
-
How to create tx2gene object with gene names from annotation
written by poecile.pal •
-
I have a list of chromosomes, positions, gene ids, and transcript ids (from a gencode GTF file) - h…
written by claudiadast •
-
Choice of Gene Length in TPM Calculations
written by dec986
-
GTF2BED does not provide meaningful output?
written by c_u
-
Comprehensive gene annotation for drosophila melanogaster
written by newscient •
-
gene and exon number from bed file
written by bioguy24
-
What's the number at the begining of some lines of transcriptome GTF file means ?
written by jack
What does the BED file look like?
My bam file looks like this:
Without
rseqcsimply looking at tags and genes can we tell whether data is strand specific or not. If so how?can you please help me with this...for my previous comment.
No we can't tell by just looking at alignments. What result did you get from running the script?
This is PairEnd Data Fraction of reads failed to determine: 0.0000
As "1+-,1-+,2++,2--" gets more signal which means fr-firststrand. So, in HISAT2 alignment step the –rna –strandness parameter will be RF.
A colleague of mine without
rseqcby looking at tags and genes told me that the data is strand specific. But I'm curious to know how he said that.This means you have a standard (dUTP-based) strand-specific library. They may have done it looking at the
XS:Atags. Every aligner does not do that so safer to check with the script.okay. you mean in the above bam it is seen that
XS:Ais given-so this tells the data is strand specific? Am I right? and how based on genes we can say that?