search

hartleys / JunctionSeq

The JunctionSeq R package is a powerful tool for testing and visualizing differential usage of exons and splice junctions in next-generation, high-throughput RNA-Seq experiments.
28 stars 16 forks source link

Annotation file appears to be missing information! Are you sure you are using the correct flattened annotation file, created by the QoRTs QC command? #35

Closed kaanokay closed 3 years ago

kaanokay commented 6 years ago

Hello, I have a problem that about number of rows of GFF file and Count file. I used hg38.gtf file that was downloaded from ENSEMBL. When I do QC.spliceJunctionAndExonCounts.forJunctionSeq file with using hg38.gtf annotation file, there isn't any problem. Then I run mergeNovelSplices command and JS.geo.size.factor file was done by this command. JS.geo.size.factor file has size factors number for each samples. Yet, when I run "runJunctionSeqAnalyses" function on R, I take a error that say: FATAL ERROR! Annotation file appears to be missing information! More information below: Count files have: 1026182 rows. Flat GFF file has: 682212 rows. WARNING: number of count file rows does not match flattned GFF file! One of the files may have been truncated, or may have been created using different strandedness rules. Missing information on 169490/502864 features. A few of the missing feature ID's: ["ENSG00000177311:E011","ENSG00000177311:E012","ENSG00000177311:E013","ENSG00000177311:E014","ENSG00000177311:E015","ENSG00000177311:E016","ENSG00000177311:E017","ENSG00000177311:E018","ENSG00000177311:E019","ENSG00000177311:E020"] Error in readJunctionSeqCounts(countfiles = as.character(sample.files), : FATAL ERROR! Annotation file appears to be missing information! Are you sure you are using the correct flattened annotation file, created by the QoRTs QC command?

Additionally, when I use hg38.gtf file that was downloaded as the UCSC format, then JS.geo.size.factor file has NaN values for each samples and I take error in R. My dataset is single ended and I think it is unstranded. ( Illumina Genome Analyzer II was used to perform RNA-seq library, Standard library preparation protocol for high throughput DNA sequencing on Illumina Genome Analyzer was performed....)

My all commands are below step by step; 1-) I performed this command for all bam files. Each samples has one bam file. ~$ java -jar Desktop/QoRTs.jar QC --singleEnded --keepMultiMapped --runFunctions writeKnownSplices,writeNovelSplices,writeSpliceExon --minMAPQ 20 Desktop/inputData/bam_files/input.bam Desktop/inputData/annoFiles/hg38_annotated_ENSEMBL.gtf Desktop/outputData/output_bam_file

2-) ~$ java -jar Desktop/QoRTs.jar makeFlatGff --verbose Desktop/inputData/annoFiles/hg38_annotated_ENSEMBL.gtf Desktop/inputData/annoFiles/JunctionSeq.flat.gff.gz

3-) ~$ java -jar Desktop/QoRTs.jar mergeNovelSplices --verbose --minCount 6 Desktop/outpuData Desktop/inputData/annoFiles/decoder.bySample.txt Desktop/inputData/annoFiles/hg38_annotated_ENSEMBL.gtf Desktop/outputData

4-) R command; JunctionSeqCountSet= runJunctionSeqAnalyses(sample.files = countFiles.noNovel, sample.names = phenotype.data$sample.ID, condition = factor(phenotype.data$group.ID), flat.gff.file = gff.file, nCores = 1, analysis.type = "junctionsAndExons")

Splice-Exon counts and GFF file appearance at below

zcat QC.spliceJunctionAndExonCounts.forJunctionSeq.txt.gz | head -n20

ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:A000 59 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E001 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E002 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E003 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E004 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E005 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E006 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E007 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E008 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E009 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E010 0 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E011 1 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E012 1 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E013 13 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E014 6 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E015 2 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E016 3 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E017 15 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E018 4 ENSG00000223972+ENSG00000227232+ENSG00000243485+ENSG00000284332:E019 6

zcat JunctionSeq.flat2.gff.gz | head -n20

chr1 ScalaUtils aggregate_gene 11874 29370 . . . gene_id NR_024540+NR_046018; tx_set NR_024540+NR_046018; num 000; aggregateGeneStrand .; geneCt 2; tx_ct 2; tx_strands -,+; tx_cds_spans 0_0,0_0 chr1 ScalaUtils exonic_part 11874 12227 . . . gene_id NR_024540+NR_046018; tx_set NR_046018; num 001; gene_set NR_046018 chr1 ScalaUtils exonic_part 12613 12721 . . . gene_id NR_024540+NR_046018; tx_set NR_046018; num 002; gene_set NR_046018 chr1 ScalaUtils exonic_part 13221 14361 . . . gene_id NR_024540+NR_046018; tx_set NR_046018; num 003; gene_set NR_046018 chr1 ScalaUtils exonic_part 14362 14409 . . . gene_id NR_024540+NR_046018; tx_set NR_024540+NR_046018; num 004; gene_set NR_024540+NR_046018 chr1 ScalaUtils exonic_part 14410 14829 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 005; gene_set NR_024540 chr1 ScalaUtils exonic_part 14970 15038 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 006; gene_set NR_024540 chr1 ScalaUtils exonic_part 15796 15947 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 007; gene_set NR_024540 chr1 ScalaUtils exonic_part 16607 16765 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 008; gene_set NR_024540 chr1 ScalaUtils exonic_part 16858 17055 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 009; gene_set NR_024540 chr1 ScalaUtils exonic_part 17233 17368 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 010; gene_set NR_024540 chr1 ScalaUtils exonic_part 17606 17742 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 011; gene_set NR_024540 chr1 ScalaUtils exonic_part 17915 18061 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 012; gene_set NR_024540 chr1 ScalaUtils exonic_part 18268 18366 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 013; gene_set NR_024540 chr1 ScalaUtils exonic_part 24738 24891 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 014; gene_set NR_024540 chr1 ScalaUtils exonic_part 29321 29370 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 015; gene_set NR_024540 chr1 ScalaUtils splice_site 12228 12612 . . . gene_id NR_024540+NR_046018; tx_set NR_046018; num 016; gene_set NR_046018 chr1 ScalaUtils splice_site 12722 13220 . . . gene_id NR_024540+NR_046018; tx_set NR_046018; num 017; gene_set NR_046018 chr1 ScalaUtils splice_site 14830 14969 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 018; gene_set NR_024540 chr1 ScalaUtils splice_site 15039 15795 . . . gene_id NR_024540+NR_046018; tx_set NR_024540; num 019; gene_set NR_024540

Can you help to me about this issue, thank a lot

hartleys commented 5 years ago

So it looks to me at first glance that "JunctionSeq.flat2.gff.gz" was created using a different (non-ensembl) original GTF. Note that the geneID's aren't ensembl ID's.

Did you intend to use "JunctionSeq.flat.gff.gz" instead?

But it looks like some of the features were matched, which is odd. Like maybe you were using different version of the ensembl annotation?

Have you tried grepping the example gene that it prints out in the error message?

 zcat JunctionSeq.flat2.gff.gz | grep ENSG00000177311

What does that give you?

kaanokay commented 5 years ago

Thanks for your answer I run this command again; ~$ java -jar Desktop/QoRTs.jar makeFlatGff --verbose Desktop/inputData/annoFiles/hg38_annotated_ENSEMBL.gtf Desktop/inputData/annoFiles/JunctionSeq.flat.gff.gz

When I run zcat JunctionSeq.flat.gff.gz | grep ENSG00000177311 command in terminal, result is ;

gzip: JunctionSeq.flat.gff.gz: unexpected end of file 3 ScalaUtils aggregate_gene 141324213 141449792 . .. gene_id ENSG00000177311; tx_set ENST00000321464+ENST00000441582+ENST00000503809+ENST00000504673+ENST00000506623+ENST00000507657+ENST00000507722+ENST00000509813+ENST00000509842+ENST00000509883+ENST00000510338+ENST00000510726+ENST00000512276+ENST00000512327+ENST00000512769+ENST00000513249+ENST00000513258+ENST00000513570+ENST00000513619+ENST00000514251+ENST00000636289+ENST00000637056+ENST00000637706; num 000; aggregateGeneStrand +; geneCt 1; tx_ct 23; tx_strands +,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+,+; tx_cds_spans 141442389_141445973,141442389_141445973,141442389_141442491,141442389_141442670,0_0,0_0,141442389_141442744,141442389_141442722,141442389_141442840,141442389_141444701,141442389_141442826,141442389_141443020,0_0,0_0,0_0,0_0,141442389_141442713,141442389_141442597,0_0,141442389_141445973,141442389_141444565,141442389_141445973,0_0 3 ScalaUtils exonic_part 141324213 141324456 . .. gene_id ENSG00000177311; tx_set ENST00000509842; num 001; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141368304 141368323 . .. gene_id ENSG00000177311; tx_set ENST00000513619; num 002; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141368324 141368443 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276; num 003; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141368444 141368489 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276+ENST00000637706; num 004; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141368490 141368620 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276+ENST00000637706+ENST00000507657; num 005; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141368621 141368804 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276+ENST00000637706+ENST00000509842+ENST00000507657; num 006; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141369872 141369946 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276+ENST00000637706+ENST00000509842+ENST00000507657; num 007; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141381425 141381487 . .. gene_id ENSG00000177311; tx_set ENST00000513619+ENST00000512276+ENST00000637706+ENST00000509842+ENST00000507657; num 008; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141383660 141383709 . .. gene_id ENSG00000177311; tx_set ENST00000637706; num 009; gene_set ENSG00000177311 3 ScalaUtils exonic_part 141384810 141384875 . .. gene_id ENSG00000177311; tx_set ENST00000507722; num 010; gene_set ENSG00000177311

İf my annotation file is not appropriate for the analysis, then can you send link which has appropriate annotation file(for human) for this analysis because I want to analyze exon usage and junction usage for my RNA-Seq data. There are many links which have ensembl annotation file, but I don't know which one is appropriate.

Also, I aligned my fastq file using Hisat2 instead of RNA-STAR. Is this problem ?

Thanks for your answer

Steve Hartley notifications@github.com, 20 Eyl 2018 Per, 19:04 tarihinde şunu yazdı:

So it looks to me at first glance that "JunctionSeq.flat2.gff.gz" was created using a different (non-ensembl) original GTF. Note that the geneID's aren't ensembl ID's.

Did you intend to use "JunctionSeq.flat.gff.gz" instead?

But it looks like some of the features were matched, which is odd. Like maybe you were using different version of the ensembl annotation?

Have you tried grepping the example gene that it prints out in the error message?

zcat JunctionSeq.flat2.gff.gz | grep ENSG00000177311

What does that give you?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHub https://github.com/hartleys/JunctionSeq/issues/35#issuecomment-423238971, or mute the thread https://github.com/notifications/unsubscribe-auth/AoR12T7Sd4tMIXSGz5j40_SyL1zHLWUyks5uc7yQgaJpZM4V6hS9 .

-- Kaan Okay Oktay Laboratory Neuro-Genomics Group Izmir International Biomedicine and Genome Center (IBG-Izmir) Dokuz Eylul University