[apredeus]

There are many possible choices for annotations: RefSeq, UCSC genes, ENSEMBL. I would recommend Gencode annotations, which are based on ENSEMBL, are very comprehensive, and are used by ENCODE.
You can simply add these annotations GTF files to the mRNA GTF files, making sure transcript IDs are distinct. This will increase the number of junctions in your database, which is usually beneficial. Again, depending on the quality of mRNA alignments you may need to filter them for junctions with very short introns, and, possibly, other artifacts. A simple way to do it is to start generating the genome with mRNA GTF file, then filter suspicious junctions from jdbList.out.tab file, and then re-generate genome feeding the filtered junctions with --sjdbFileChrStartEnd (you can include the annotations GTF file at the same time).

Yes, and in order to avoid Gencode/UCSC scaffold naming differences, I will just use the Gencode GTF with no "random" annotations. Sounds great. Thanks again for the answers!

[apredeus]

I recommend including the *GL* or *gl* marked "unplaced" scaffolds. There are some rRNAs on these scaffolds from which large number of reads may originate, especially if the ribo-depletion did not work well.

On the other hand, *hap* scaffolds represent haplotypes. Some reads will map equally well to multiple haplotypes, and thus will be marked as "multi-mappers", which is not a desired behavior in most cases. Unless you need this kind of halpotype-aware mapping, I do not recommend including them.

Yes, very logical. Funny, after I posted this question, this whole story came up:

http://seqanswers.com/forums/showthread.php?p=136425

So, definitely including the extra scaffolds!

[alexdobin]

Also, I meant to ask - do you normally include random chromosomes and "hap" chromosomes from human genome into the overall index? I know it shouldn't make a huge difference, however there's quite a few transcripts that are mapped to these.

Thank you!

I recommend including the *GL* or *gl* marked "unplaced" scaffolds. There are some rRNAs on these scaffolds from which large number of reads may originate, especially if the ribo-depletion did not work well.

On the other hand, *hap* scaffolds represent haplotypes. Some reads will map equally well to multiple haplotypes, and thus will be marked as "multi-mappers", which is not a desired behavior in most cases. Unless you need this kind of halpotype-aware mapping, I do not recommend including them.

[alexdobin]

Great, thank you very much for the informative answer. I'll make sure to filter out the ultra-short introns.



What would "standard" annotation be in this case? RefSeq? It's just that we have always used mRNA collection for both humans and mice (there's about 1.5 mil for mm9 and 2.5 or so for hs19), I'm not sure what else do people use.

There are many possible choices for annotations: RefSeq, UCSC genes, ENSEMBL. I would recommend Gencode annotations, which are based on ENSEMBL, are very comprehensive, and are used by ENCODE.
You can simply add these annotations GTF files to the mRNA GTF files, making sure transcript IDs are distinct. This will increase the number of junctions in your database, which is usually beneficial. Again, depending on the quality of mRNA alignments you may need to filter them for junctions with very short introns, and, possibly, other artifacts. A simple way to do it is to start generating the genome with mRNA GTF file, then filter suspicious junctions from jdbList.out.tab file, and then re-generate genome feeding the filtered junctions with --sjdbFileChrStartEnd (you can include the annotations GTF file at the same time).

[apredeus]

Also, I meant to ask - do you normally include random chromosomes and "hap" chromosomes from human genome into the overall index? I know it shouldn't make a huge difference, however there's quite a few transcripts that are mapped to these.

Thank you!

[apredeus]

Great, thank you very much for the informative answer. I'll make sure to filter out the ultra-short introns.

If you want to use mRNA GTF files instead or in addition to standard annotations, I would recommend checking the splice junctions in this file for very short introns, and filtering them out - please see this post.

What would "standard" annotation be in this case? RefSeq? It's just that we have always used mRNA collection for both humans and mice (there's about 1.5 mil for mm9 and 2.5 or so for hs19), I'm not sure what else do people use.

[alexdobin]

Alex, thank you for the great tool - STAR is indeed very impressive!

Would you post a command that one would use to generate a mixed genome for mouse and human? We normally align to mouse and human genomes using GTF files of annotated mRNAs. Can one incorporate both of those GTF files into a combined mouse+human index?

Hi Alex,

to generate combined mouse/human genome with STAR you would need to modify slightly your fasta and gtf files:
1. Modify chromosomes names so that mouse and human chromosomes have distinct names, e.g. chr1h/chr1m etc. In the FASTA files you need to make these modifications in all sequences name lines (i.e. starting with ">"). In GTF files you would need to modify all chromosome names in field 1.
2. Make sure that the transcript_id in GTF files are distinct for mouse and human. This is usually the case, for instance, Gencode has "ENSMUSTxxxxx" for mouse and "ENSTxxxxx" for human.
3. Concatenate GTF files for mouse and human into a single GTF file
4. Run genome generation with
STAR --runMode genomeGenerate --runThreadN 12 --genomeDir ./ --genomeFastaFiles /path/to/human.fa /path/to/mouse.fa --sjdbGTFfile /path/to/mouse_human.gtf --sjdbOverhang 100

If you want to use mRNA GTF files instead or in addition to standard annotations, I would recommend checking the splice junctions in this file for very short introns, and filtering them out - please see this post.

Cheers
Alex

[apredeus]

We used STAR in two mixed-species settings: human + viruses (~4,500 viruses from NCBI) and human + mouse, and in both cases it worked well as far as we could tell. You have to keep an eye on RAM, since STAR would need ~10*TotalGenomeSize bytes of RAM (~50GB for human+mouse), and if you have a large number of small chromosomes/scaffolds/contigs, you would need to reduce --genomeChrBinNbits as I explained in the previous post

Alex, thank you for the great tool - STAR is indeed very impressive!

Would you post a command that one would use to generate a mixed genome for mouse and human? We normally align to mouse and human genomes using GTF files of annotated mRNAs. Can one incorporate both of those GTF files into a combined mouse+human index?

[dpryan]

Hey Devon,

Its turns it is not difficult since a group of individual from Case Western Reserve University, Cleveland, OH published a paper on a program they developed called LoQuM which does exactly what I wanted. I have not tried the program yet but here is the title of article if you would like to read if yourself

"Accurate estimation of short read mapping quality for next-generation genome sequencing"

Thanks,
Nino

Interesting, I'll have to give that paper a read, thanks!

[wildtypegoose]

STAR will use a 1-2 GB of memory per process for temporary storage and I/O buffers, however, the ~25GB of genome files are shared. The "used" memory reported by top includes "cached", and it's hard to determine how much physical RAM the process is actually using.

Hi Alex,
You are right: I was able to run 7 star jobs simultaneously on our server with 128GB of memory. I used the LoadAndExit option of genomeLoad flag to first load the genome in shared memory, and then used the LoadAndRemove for all the simultaneous STAR jobs.

Although the processes went well on memory usage, but I noticed that server at times became irresponsible due to a lot of I/O (as shown by D state for many of the STAR processes in "top" output). Any suggestion/s to avoid this bottleneck?

Thanks a lot!

[Nino]

Hey Devon,

Its turns it is not difficult since a group of individual from Case Western Reserve University, Cleveland, OH published a paper on a program they developed called LoQuM which does exactly what I wanted. I have not tried the program yet but here is the title of article if you would like to read if yourself

"Accurate estimation of short read mapping quality for next-generation genome sequencing"

Thanks,
Nino

[dpryan]

Hello

I am happy with STAR but not too happy with the MQ scores it is difficult to filter out reads based on their mapping quality score when you only have 255,3,2,1,0 (I believe 0 is one of them I forget). Does anyone know of any program which can convert these 5 or 4 values based on the CIGAR information of the read to phred scale values?

Thanks,
Nino

That turns out to be a surprisingly difficult thing to do, as you often end up needing to realign everything so that you know how many second-best alignments there are and what their score is (unless one of STAR's more verbose output modes provides this).

[Nino]

Hello

I am happy with STAR but not too happy with the MQ scores it is difficult to filter out reads based on their mapping quality score when you only have 255,3,2,1,0 (I believe 0 is one of them I forget). Does anyone know of any program which can convert these 5 or 4 values based on the CIGAR information of the read to phred scale values?

Thanks,
Nino

[alexdobin]

The global "used" memory (as reported by top) increases linearly if I run another STAR process, which made me think that probably the shared mem option is not working correctly. I'll try your suggestion of running >5 jobs once the server is free.

Thanks for your input!

STAR will use a 1-2 GB of memory per process for temporary storage and I/O buffers, however, the ~25GB of genome files are shared. The "used" memory reported by top includes "cached", and it's hard to determine how much physical RAM the process is actually using.

[wildtypegoose]

It seems to me that the shared memory is working fine - since ipcs shows two instances attached to the same shared memory piece. What you see in 'top' is memory usage per process - and 25GB out of 26.6GB are shared for each process. To test it you can try to run >5 jobs at the same time - they would not be able to run without sharing memory.

The global "used" memory (as reported by top) increases linearly if I run another STAR process, which made me think that probably the shared mem option is not working correctly. I'll try your suggestion of running >5 jobs once the server is free.

Thanks for your input!