[thomasblomquist]

If I understand it correctly, you need to match short sub-sequences of your long reads to a set of reference short sequences. For starters I would suggest the following parameters:
--outFilterMismatchNmax N (N=number of mismatches you wish to tolerate)
--outFilterScoreMinOverLread 0
--outFilterMatchNminOverLread 0
--outFilterMatchNmin L (L=shortest of the reference sequences)

Please let me know whether it works for you.

In essence. Since I generally know where the alleles are, the sequence reads are trimmed to include the variant loci and ~15 bases of sequencing read on either side of the expected loci(s). So, I'm matching ~40 base sequences to a reference library of chromosomes that are ~20-30 bases.

Thanks again for your help. I will try the modifications later tonight. -Tom

[alexdobin]

Thank you. I do tend to go about things in a non-standard way . I am performing targeted resequencing of an RNA-seq library. Also known as amplicon sequencing. It is also semi-quantitative. My end goal being, I want to count the number of times a given sequence (e.g. for a gene target), or genetic variation (alleles from 1 base substitution up to 6 bases substituted). I want to be able to assign each sequence read the "name" of the consensus sequence it best matches.

So, my reference genome is usually anywhere from 10-1000 short sequences of ~10-40 bases in length.

Thank you for your time. I look forward to seeing how your software works in my studies.

Regards,

-Tom Blomquist

If I understand it correctly, you need to match short sub-sequences of your long reads to a set of reference short sequences. For starters I would suggest the following parameters:
--outFilterMismatchNmax N (N=number of mismatches you wish to tolerate)
--outFilterScoreMinOverLread 0
--outFilterMatchNminOverLread 0
--outFilterMatchNmin L (L=shortest of the reference sequences)

Please let me know whether it works for you.

[thomasblomquist]

Tom,
please use the following parameters:
--genomeChrBinNbits 6 --genomeSAindexNbases 4
I was able to run successful genome generation step for your small reference.
It looks like you want to use STAR in a quite non-standard way.
If you explain it a bit more I could suggest the parameters for the mapping step.

Thank you. I do tend to go about things in a non-standard way . I am performing targeted resequencing of an RNA-seq library. Also known as amplicon sequencing. It is also semi-quantitative. My end goal being, I want to count the number of times a given sequence (e.g. for a gene target), or genetic variation (alleles from 1 base substitution up to 6 bases substituted). I want to be able to assign each sequence read the "name" of the consensus sequence it best matches.

So, my reference genome is usually anywhere from 10-1000 short sequences of ~10-40 bases in length.

Thank you for your time. I look forward to seeing how your software works in my studies.

Regards,

-Tom Blomquist

[alexdobin]

small reference

I'm having trouble generating a small reference genome for running STAR. Below is the command, Error output and the input "genome" file. Your thoughts on how I can rectify are welcome. At present, I am using BFAST and it works reasonably well, but am looking for something a little easier for automatically generating seeds for parsing sequencing reads into bins. Thank you for your help and time.

Regards,

-Tom Blomquist

University of Toledo, Ohio

Tom,
please use the following parameters:
--genomeChrBinNbits 6 --genomeSAindexNbases 4
I was able to run successful genome generation step for your small reference.
It looks like you want to use STAR in a quite non-standard way.
If you explain it a bit more I could suggest the parameters for the mapping step.

[thomasblomquist]

Trouble generating a small genome suffix array index...

I'm having trouble generating a small reference genome for running STAR. Below is the command, Error output and the input "genome" file. Your thoughts on how I can rectify are welcome. At present, I am using BFAST and it works reasonably well, but am looking for something a little easier for automatically generating seeds for parsing sequencing reads into bins. Thank you for your help and time.

Regards,

-Tom Blomquist

University of Toledo, Ohio



Command:

./STAR --runMode genomeGenerate --genomeDir ./ --genomeFastaFiles ./barcode.fa --runThreadN 8



Error Output:

Feb 19 09:39:34 ..... Started STAR run
Feb 19 09:39:34 ... Starting to generate Genome files
Feb 19 09:39:34 ... starting to sort Suffix Array. This may take a long time...
Feb 19 09:39:34 ... sorting Suffix Array chunks and saving them to disk...
Feb 19 09:39:34 ... loading chunks from disk, packing SA...
Feb 19 09:39:34 ... writing Suffix Array to disk ...
Feb 19 09:39:34 ... Finished generating suffix array
Feb 19 09:39:34 ... starting to generate Suffix Array index...

BUG: next index is smaller than previous, EXITING



./barcode.fa is the following -->

>ATGC
TTTTCATGCGATCAGGCGTCTGTCGTGCTC

>CAGT
TTTTCCAGTGATCAGGCGTCTGTCGTGCTC

>GACA
TTTTCGACAGATCAGGCGTCTGTCGTGCTC

>TGTG
TTTTCTGTGGATCAGGCGTCTGTCGTGCTC

>ACAC
TTTTCACACGATCAGGCGTCTGTCGTGCTC

>TACG
TTTTCTACGGATCAGGCGTCTGTCGTGCTC

>GCTA
TTTTCGCTAGATCAGGCGTCTGTCGTGCTC

>CATA
TTTTCCATAGATCAGGCGTCTGTCGTGCTC

>TAGA
TTTTCTAGAGATCAGGCGTCTGTCGTGCTC

>GACT
TTTTCGACTGATCAGGCGTCTGTCGTGCTC

>CTGA
TTTTCCTGAGATCAGGCGTCTGTCGTGCTC

>TGCT
TTTTCTGCTGATCAGGCGTCTGTCGTGCTC

[alexdobin]

large number of contigs

I've also had success in working around this problem by creating dummy scaffolds. I cat all the contigs into a single big fasta entry with a couple hundred 'X's separating each contig to (hopefully) prevent spurious alignments from spanning contigs. This made a huge difference in the amount of memory required.

This is a great idea. It will solve both the RAM problem and the slowdown problem I was talking about. However, you would need to post-process your alignments to extract alignments' coordinates within the real contigs from the coordinates within "super-contigs". Also, while Xs or Ns separating the contigs will prevent STAR from "extending" the alignment into adjacent contigs, it will not prevent splicing between contigs. So at the post-processing step you would also need to filter out alignments that span more than one real contig.

[cram]

As @NGSfan pointed out, this is indeed RAM problem. STAR bins genome sequence in a way that each chromosome (contig) starts at a new bin, which creates an overhead of Nchromosomes*BinSize, where BinSize=2^genomeChrBinNbits. By default, --genomeChrBinNbits = 18,
so BinSize=2^18~256kb, so with 300,000 contigs you would need ~75GB of RAM - that's what likely killed your job.

I suggest that you try a much smaller value of --genomeChrBinNbits 12. This would require just a few GB of RAM and should allow you to generate the genome files. I have not tried STAR with more than 50,000 contigs, and I suspect there might be significant slowdown in the mapping speed when the number of contigs is too big.

I've also had success in working around this problem by creating dummy scaffolds. I cat all the contigs into a single big fasta entry with a couple hundred 'X's separating each contig to (hopefully) prevent spurious alignments from spanning contigs. This made a huge difference in the amount of memory required.

[alexdobin]

mixed datasets

This looks great, I think I'll try it out! Have you tested STAR on mixed RNA datasets? i.e. RNA-Seq on a sample containing both human and bacterial/viral RNA? Thats what I'm working with at the moment and Tophat just isn't cutting it, there are millions of reads that are identified as being human and bacterial (depending on which aligning step I run first)

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.

For mixed-species mapping I strongly agree with other users who recommended mapping to a combined genome rather than a 2-step mapping. When a mapper aligns a read to combined genome it considers various possible alignments and picks the best one. This reduces the false positives and negatives for unique calls to each species. On the other hand, in the 2-step method you force alignments to one of the species in the 1st step, so you results will be strongly biased towards the 1st species.

[alexdobin]

large number of contigs

I have problems generating a new genome.
My genome is more than 300.000 contigs and ~300MB in size.

Feb 18 11:25:11 ..... Started STAR run
Feb 18 11:25:11 ... Starting to generate Genome files
/var/spool/slurmd/job1306271/slurm_script: line 9: 12129 Killed

As @NGSfan pointed out, this is indeed RAM problem. STAR bins genome sequence in a way that each chromosome (contig) starts at a new bin, which creates an overhead of Nchromosomes*BinSize, where BinSize=2^genomeChrBinNbits. By default, --genomeChrBinNbits = 18,
so BinSize=2^18~256kb, so with 300,000 contigs you would need ~75GB of RAM - that's what likely killed your job.

I suggest that you try a much smaller value of --genomeChrBinNbits 12. This would require just a few GB of RAM and should allow you to generate the genome files. I have not tried STAR with more than 50,000 contigs, and I suspect there might be significant slowdown in the mapping speed when the number of contigs is too big.

[alexdobin]

read count

@iris_aurelia:
My plan is to count the reads per gene within STAR algortihm as they are mapped, so it will be as multi-threaded as STAR itself.
Multi-threading on a per-chromosome basis is a good idea. I am going to implement a multi-threaded BAM sorting using this kind of parallelization.

@NGSfan:
Using Cufflinks "coverage" for calculating gene counts is a nice trick. We might as well do it since we are typically running both Cufflinks and HT-seq on our datasets. In our experience Cufflinks' speed is acceptable for A+ datasets, however, it slows down considerably and often gets stuck on A- and Total RNA samples.

[bob-loblaw]

This looks great, I think I'll try it out! Have you tested STAR on mixed RNA datasets? i.e. RNA-Seq on a sample containing both human and bacterial/viral RNA? Thats what I'm working with at the moment and Tophat just isn't cutting it, there are millions of reads that are identified as being human and bacterial (depending on which aligning step I run first)

[JonB]

NGSfan:
Thanks. I think I had 60 GB of memory last time. I will try some more

[NGSfan]

In my experience Cufflinks (using the multithreaded option) is even slower than HTseq-count.
We are doing some experimenting creating our own multithreaded count script which outputs exact the same data as HT-seq count. In order to do this multithreaded we are running each chromosome separately, which might be an idea for STAR?

That's good to know - I have been using 24 threads which is pretty fast, but I did try it on 8 threads and noticed it is not much faster at that level, especially if you include options for multi-read correction and fragment bias correction.

I would support any attempt to multi-thread HTseq-count, let us know when you have a script. You can share it on code.google.com, it's where I put my most useful scripts.

[NGSfan]

I have problems generating a new genome.
My genome is more than 300.000 contigs and ~300MB in size.

Feb 18 11:25:11 ..... Started STAR run
Feb 18 11:25:11 ... Starting to generate Genome files
/var/spool/slurmd/job1306271/slurm_script: line 9: 12129 Killed

When the genome has a lot of contigs, I noticed that STAR needs a lot more memory... how much do you have? I would recommend finding a machine with 256gb RAM if possible.

[iris_aurelia]

That is a great feature, Alex. I am looking forward to it.

I currently convert cufflinks output into read counts for DESeq.

You read in the length and coverage columns from the isoforms.fpkm_tracking file as well as the read length of your reads and then caclulate:

my $reads = $length * $coverage / $readlength;

I have compared the values to htseq-count and found them to be nearly the same.

The advantage over htseq is you can multi-thread cufflinks.

In my experience Cufflinks (using the multithreaded option) is even slower than HTseq-count.
We are doing some experimenting creating our own multithreaded count script which outputs exact the same data as HT-seq count. In order to do this multithreaded we are running each chromosome separately, which might be an idea for STAR?