[serash]

Hi,
I'm currently doing 2 pass alignment with the --twopass1readsN option, which adds the 1st pass junctions into genome in memory. I'm curious to know if this works when a genome in memory i used (--genomeLoad LoadAndKeep) from previous runs. And will this give the correct results if I run additional 2 pass alignments with this genome in memory or do I need to reload the genome for each run?

Thanks!

[coryfunk]

My mistake. Thanks!

[Brian Bushnell]

0x4 means 'unmapped'. If no flags are set, that means it's a single-ended read mapped to the plus strand.

[dpryan]

No, not setting that bit means the read is mapped.
Edit:I'll add that samtools view -c -F 4 foo.bam will count the mapped reads. Using -f 4 will count the unmapped reads, which will be 0.

[coryfunk]

That's a bit of an issue because the default of not having a flag set is to call it "unmapped".

Running:

samtools view -c -F 4 <file.name>

counts all the reads as unmapped. Downstream tools that rely on that flag are thus rendered incompatible.

[dpryan]

STAR doesn't output unmapped reads in the SAM file (at least not by default), so bit 4 shouldn't be set there for any record.

[coryfunk]

In using STAR for my alignments, I'm subsequently using samtools for some basic evaluations (and GATK after that to generate a vcf file). An issue I'm having is that even though STAR is correctly including the mapping location for each read in the sam file, it is not setting the flag that says if a read is mapped or not (-f 4). As such, when I use samtools to count the number of mapped reads, it's always zero. Samtools accurately counts the number of total reads and unmapped reads.

Is this an issue with STAR or possibly something I'm doing incorrectly?

[ndaniel]

Hi Emma,

these long-gap splices, often connecting adjacent genes, are somewhat common in RNA-seq data. It's hard to say whether they are biochemically real "read-through transcription" events, or some kind of wet-lab or mapping artifacts.
They would be clearly mapping artifacts if "better" alignments of these sequences can be found, however, BLATing or BLASTing them did not result in any better alignments.
One way to get rid of them is to completely prohibit long gaps with --alignIntronMax N, which would prohibit any gap longer than N (by default this is ~600000). However, if you make this too small, say 100000, you may miss a number of valid junctions, as mammalian introns can be hundred of kilobases long.
A better approach is filter out long-gap alignments supported by too few reads, e.g. :
--outFilterType BySJout --outSJfilterIntronMaxVsReadN 10000 20000 50000 100000
This would only allow unannotated junctions <=10kb supported by >=1 spliced read, <=20kb supported by >=2 reads, <=50kb by >= 3 reads, <=100kb by >=4 reads.

This looks like a very bad idea. There are plenty of fusion genes which are known to be "adjacent" but still are genuine fusion genes.

For example, the famous fusion gene FGFR3-TACC3 (for more info about it please Google it) has the genes FGFR3 and TACC3 as adjacent and the distance between them is less than 50kb. If this fusion is expressed moderately there could be, for example, 2 reads supporting the fusion and it will be discarded according to the above criteria.

Out there are several fusion finders which can handle this kind of cases with very low false positives rate, e.g. https://code.google.com/p/fusioncatcher/wiki/comparison

[alexdobin]

Hi,
Trying to generate genome from Rice reference and I get the following error,have tried several STAR patches available:

biostat1@biostat[STAR_2.3.1z10] ./STAR --runMode genomeGenerate --genomeDir IRGSP_genome --genomeFastaFiles L1_1.fq L1_2.fq
Jul 17 10:55:11 ..... Started STAR run
Jul 17 10:55:11 ... Starting to generate Genome files
terminate called after throwing an instance of 'std:ut_of_range'
what(): vector::_M_range_check
zsh: abort ./STAR --runMode genomeGenerate --genomeDir IRGSP_genome --genomeFastaFiles

Any ideas,
Thanks!

Hi @kjusto,
Is it possible that your L1_1.fq L1_2.fq files are not fasta files - they seem to have FASTQ extension?
If they are valid fasta, please send me a link to them, and I will try to figure out the problem.
Cheers
Alex

[alexdobin]

Hi Brian,
No, I didn't use a gtf to do the mapping in this case. Plus the mapped read corresponds to a known intron (of a short gene) on one side but a random intergenic region way after the end of the gene of which it starts in ( at least according to UCSC) on the other side. And the 200kb overlaps with 4 other known genes too. So to my eyes thats definitely a mapping error too. The question now is how to filter those out (because they are quite a few of them).

Hi Emma,

these long-gap splices, often connecting adjacent genes, are somewhat common in RNA-seq data. It's hard to say whether they are biochemically real "read-through transcription" events, or some kind of wet-lab or mapping artifacts.
They would be clearly mapping artifacts if "better" alignments of these sequences can be found, however, BLATing or BLASTing them did not result in any better alignments.
One way to get rid of them is to completely prohibit long gaps with --alignIntronMax N, which would prohibit any gap longer than N (by default this is ~600000). However, if you make this too small, say 100000, you may miss a number of valid junctions, as mammalian introns can be hundred of kilobases long.
A better approach is filter out long-gap alignments supported by too few reads, e.g. :
--outFilterType BySJout --outSJfilterIntronMaxVsReadN 10000 20000 50000 100000
This would only allow unannotated junctions <=10kb supported by >=1 spliced read, <=20kb supported by >=2 reads, <=50kb by >= 3 reads, <=10kb by >=4 reads.

There is more discussion on this type of filtering in this post.

Cheers
Alex

[kjusto]

Hi,
Trying to generate genome from Rice reference and I get the following error,have tried several STAR patches available:

biostat1@biostat[STAR_2.3.1z10] ./STAR --runMode genomeGenerate --genomeDir IRGSP_genome --genomeFastaFiles L1_1.fq L1_2.fq
Jul 17 10:55:11 ..... Started STAR run
Jul 17 10:55:11 ... Starting to generate Genome files
terminate called after throwing an instance of 'std:ut_of_range'
what(): vector::_M_range_check
zsh: abort ./STAR --runMode genomeGenerate --genomeDir IRGSP_genome --genomeFastaFiles

Any ideas,
Thanks!

[emmanouela]

However, if that 200kbp corresponds exactly to a known intron in the GTF file, and only occurs when using the GTF file, it's probably OK. Does it?

Hi Brian,
No, I didn't use a gtf to do the mapping in this case. Plus the mapped read corresponds to a known intron (of a short gene) on one side but a random intergenic region way after the end of the gene of which it starts in ( at least according to UCSC) on the other side. And the 200kb overlaps with 4 other known genes too. So to my eyes thats definitely a mapping error too. The question now is how to filter those out (because they are quite a few of them).

[GenoMax]

Don't think Alex provides 32-bit binaries. If you have a large genome (~ human) 32-bit may not work.

Build from source if you must have 32-bit: https://code.google.com/p/rna-star/d...e.tgz&can=2&q=

[kjusto]

Pre-compiled linux binary is available here: https://code.google.com/p/rna-star/d...4.tgz&can=2&q=

Thanks for the link... got use proxies to get it though....google issues here....my question was about 32 bit linux OS,any binaries for it.

[Brian Bushnell]

Not sure about the second one, but the first one with a 200kbp deletion anchored by a 7bp of read sequence looks like a probable mapping error to me, considering that a 7bp exact match would be expected purely by chance within about 16kbp of any random location. However, if that 200kbp corresponds exactly to a known intron in the GTF file, and only occurs when using the GTF file, it's probably OK. Does it?