[alexdobin]

Hey,

I just switched from BWA to STAR for my RNAseq analysis yesterday. I got stuck with this: Mar 27 08:14:36 ..... Started STAR run
Mar 27 08:24:33 ..... Started mapping
Killed: 9

Then STAR stopped.
I have 20 of 50bp SE samples (~40M reads/sample) and the analysis was done with my poor Macmini (Core i7 / 16GB RAM).
I used Ensembl mouse genome (Mus_musculus.GRCm38.dna.primary_assembly.fa.gz and
Mus_musculus.GRCm38.79.gtf.gz) as reference.

The code I used to run the mapping:
> STAR --genomeDir ./GenomeDir --readFilesIn ./BGI_RNAseq_data_2015/01.fq,./BGI_RNAseq_data_2015/02.fq,./BGI_RNAseq_data_2015/03.fq,./BGI_RNAseq_data_2015/04.fq,./BGI_RNAseq_data_2015/05.fq,./BGI_RNAseq_data_2015/06.fq,./BGI_RNAseq_data_2015/07.fq,./BGI_RNAseq_data_2015/08.fq,./BGI_RNAseq_data_2015/09.fq,./BGI_RNAseq_data_2015/10.fq,./BGI_RNAseq_data_2015/11.fq,./BGI_RNAseq_data_2015/12.fq,./BGI_RNAseq_data_2015/13.fq,./BGI_RNAseq_data_2015/14.fq,./BGI_RNAseq_data_2015/15.fq,./BGI_RNAseq_data_2015/16.fq,./BGI_RNAseq_data_2015/17.fq,./BGI_RNAseq_data_2015/18.fq,./BGI_RNAseq_data_2015/19.fq,./BGI_RNAseq_data_2015/20.fq --runThreadN 8

I knew that STAR demands much more computing power (and memory) than BWA but I am not sure if my code is correct. I had no issues with BWA on the same Macmini though.
The log files are attached.

Hi,
to fit the mouse genome into the 16GB of RAM, you need to use at the genome generation step:
--genomeSAsparseD 2 --genomeSAindexNbases 13
Hopefully this will solve the problem.
Cheers
Alex

[neokao]

Hey,

I just switched from BWA to STAR for my RNAseq analysis yesterday. I got stuck with this: Mar 27 08:14:36 ..... Started STAR run
Mar 27 08:24:33 ..... Started mapping
Killed: 9

Then STAR stopped.
I have 20 of 50bp SE samples (~40M reads/sample) and the analysis was done with my poor Macmini (Core i7 / 16GB RAM).
I used Ensembl mouse genome (Mus_musculus.GRCm38.dna.primary_assembly.fa.gz and
Mus_musculus.GRCm38.79.gtf.gz) as reference.

The code I used to run the mapping:
> STAR --genomeDir ./GenomeDir --readFilesIn ./BGI_RNAseq_data_2015/01.fq,./BGI_RNAseq_data_2015/02.fq,./BGI_RNAseq_data_2015/03.fq,./BGI_RNAseq_data_2015/04.fq,./BGI_RNAseq_data_2015/05.fq,./BGI_RNAseq_data_2015/06.fq,./BGI_RNAseq_data_2015/07.fq,./BGI_RNAseq_data_2015/08.fq,./BGI_RNAseq_data_2015/09.fq,./BGI_RNAseq_data_2015/10.fq,./BGI_RNAseq_data_2015/11.fq,./BGI_RNAseq_data_2015/12.fq,./BGI_RNAseq_data_2015/13.fq,./BGI_RNAseq_data_2015/14.fq,./BGI_RNAseq_data_2015/15.fq,./BGI_RNAseq_data_2015/16.fq,./BGI_RNAseq_data_2015/17.fq,./BGI_RNAseq_data_2015/18.fq,./BGI_RNAseq_data_2015/19.fq,./BGI_RNAseq_data_2015/20.fq --runThreadN 8

I knew that STAR demands much more computing power (and memory) than BWA but I am not sure if my code is correct. I had no issues with BWA on the same Macmini though.
The log files are attached.

[alexdobin]

Hi,

for a 6Gb genome, you would need ~60GB of RAM, and you would need to specify --limitGenomeGenerateRAM 60000000000 (i.e. the number of bytes).

Cheers
Alex

Hi, I'm trying to use genomegenerate but every time I try to run it I have this problem:

terminate called after throwing an instance of 'std::bad_alloc'
Aborted (core dumped)

I try to use this options as well:
--runThreadN 12 --limitGenomeGenerateRAM 40

but I still get the same error.

I'm working with a genome of 6G. could that be the problem.

I really appreciate any help,

V

[vromanr_2015]

memory problems while running STAR

Hi, I'm trying to use genomegenerate but every time I try to run it I have this problem:

terminate called after throwing an instance of 'std::bad_alloc'
Aborted (core dumped)

I try to use this options as well:
--runThreadN 12 --limitGenomeGenerateRAM 40

but I still get the same error.

I'm working with a genome of 6G. could that be the problem.

I really appreciate any help,

V

[ruggedtextile]

I have narrowed the problem down to using a genome index without pre-existing junctions annotated. If I use an otherwise identical index with junctions it runs fine (takes a couple of minutes max to insert the 1st pass junctions).

I.e:

[15:49 guttea@Maul MODY] > STAR --genomeDir /mnt/nas_omics/shared_data/iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/StarAnnotatedIndex --readFilesIn data/fastq/D0-1_R1.fastq.gz data/fastq/D0-1_R2.fastq.gz --runThreadN 12 --genomeLoad NoSharedMemory --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --twopass1readsN -1 --sjdbOverhang 149
Jan 30 15:58:43 ..... Started STAR run
Jan 30 15:59:32 ..... Started 1st pass mapping
Jan 30 16:00:35 ..... Finished 1st pass mapping
Jan 30 16:02:23 ..... Finished inserting 1st pass junctions into genome
Jan 30 16:02:25 ..... Started mapping
Jan 30 16:03:43 ..... Started sorting BAM
Jan 30 16:04:04 ..... Finished successfully
[16:04 guttea@Maul MODY] > STAR --genomeDir /mnt/nas_omics/shared_data/iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/StarIndex --readFilesIn data/fastq/D0-1_R1.fastq.gz data/fastq/D0-1_R2.fastq.gz --runThreadN 12 --genomeLoad NoSharedMemory --readFilesCommand zcat --outSAMtype BAM SortedByCoordinate --twopass1readsN -1 --sjdbOverhang 149 --outFileNamePrefix NoExisitingAnnotation
Jan 30 16:06:00 ..... Started STAR run
Jan 30 16:06:15 ..... Started 1st pass mapping
Jan 30 16:07:19 ..... Finished 1st pass mapping

With this being copied at 16:15 (i.e. it has already taken ~8 minutes and last time it hung here for 18 hours using 100% CPU). I've attached the two log.out files. You can see in the NoExistingAnnotation.log.out we never get confirmation of:

Jan 30 16:01:18   Finished sorting SA indicesL nInd=50917350

Hope that helps.

[alexdobin]

I am trying to use the new (to me) twopass1readsN option in STAR, but am finding it very slow (to the point of I am not sure if STAR is crashing or not).

Using some small test data (1M PE reads) and a very small (1000) twopass1reads N parameter it takes about one minute to insert the 1st pass junctions:

STAR --genomeDir /iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/StarIndex --readFilesIn R1.fastq.gz R2.fastq.gz --runThreadN 12 --genomeLoad NoSharedMemory --outFilterMultimapNmax 1 --outSAMtype BAM SortedByCoordinate --twopass1readsN 1000 --sjdbOverhang 99
Jan 27 12:56:01 ..... Started STAR run
Jan 27 12:56:13 ..... Started 1st pass mapping
Jan 27 12:56:15 ..... Finished 1st pass mapping
Jan 27 12:57:25 ..... Finished inserting 1st pass junctions into genome
Jan 27 12:57:25 ..... Started mapping
Jan 27 12:57:41 ..... Started sorting BAM

Running on a real dataset (~40M PE reads with twopass1readsN set to -1) took >18h before I killed the job - assuming it had crashed. It was using 100% CPU that entire time however, so it may well have simply been doing the right thing - just very slowly.

Is this timing normal for this option? Manually doing the 1st pass, genome regeneration and 2nd pass would be much quicker I think. At this speed I won't be able to use the auto 2 pass mode in my normal pipeline.

Hi,

18 hours is definitely too long, 2nd pass should normally take just a bit longer than the 1st pass. Could you please send me the Log.out file of the failed run?

Cheers
Alex

[ruggedtextile]

I am trying to use the new (to me) twopass1readsN option in STAR, but am finding it very slow (to the point of I am not sure if STAR is crashing or not).

Using some small test data (1M PE reads) and a very small (1000) twopass1reads N parameter it takes about one minute to insert the 1st pass junctions:

STAR --genomeDir /iGenomes/Homo_sapiens/Ensembl/GRCh37/Sequence/StarIndex --readFilesIn R1.fastq.gz R2.fastq.gz --runThreadN 12 --genomeLoad NoSharedMemory --outFilterMultimapNmax 1 --outSAMtype BAM SortedByCoordinate --twopass1readsN 1000 --sjdbOverhang 99
Jan 27 12:56:01 ..... Started STAR run
Jan 27 12:56:13 ..... Started 1st pass mapping
Jan 27 12:56:15 ..... Finished 1st pass mapping
Jan 27 12:57:25 ..... Finished inserting 1st pass junctions into genome
Jan 27 12:57:25 ..... Started mapping
Jan 27 12:57:41 ..... Started sorting BAM

Running on a real dataset (~40M PE reads with twopass1readsN set to -1) took >18h before I killed the job - assuming it had crashed. It was using 100% CPU that entire time however, so it may well have simply been doing the right thing - just very slowly.

Is this timing normal for this option? Manually doing the 1st pass, genome regeneration and 2nd pass would be much quicker I think. At this speed I won't be able to use the auto 2 pass mode in my normal pipeline.

[alexdobin]

Hi,
Would it be wise to generate the genome with GTF file and re-generate it later with the sample specific SJ.out.tab?

I have more than 100 samples for variant analysis. Another strategy that I thing could be useful is to pull 20-30 samples together (mapping being extremely fast, it was over within 40mins), so that the discovery of novel junctions will get saturated, and then use this final SJ.out.tab to re-generate the genome, and perform sample wise alignment. This has reported ~2,00,000 total junctions, which is quite closer to the number of exons in human genome ~1,80,000

Cheers,
Amol

Hi Amol,

you are talking about what we now call a "2-pass" approach, where junctions detected in the first pass are used as annotations in the 2nd pass. This is generally the most accurate way to map the RNA-seq data. As you pointed out, there are two options:

1). Using the 1st pass in the sample specific way.
This can be done with --twopass1readsN <Nreads> option, where Nreads is the number of reads to map in the 1st pass. I recommend to use a very big number (or -1) which will map all the reads in the 1st pass (of course, in any case all the reads are re-mapped in the 2nd pass). At the moment this approach cannot be used with annotations.

2). You can collect all the junctions from the 1st mapping pass over multiple samples, and then re-generate the genome with a full collection of junctions, and re-map all reads with the new genome. This has to be done "manually" and it's a bit more cumbersome, but in the end it will give you the best sensitivity. Some filtering of the 1st pass novel junctions may be required, e.g. removing junctions on mitochondrion genomes. Some users also reported a speed decrease in the 2nd pass with millions of junctions. You will also see increase in multi-mappers, since the reads will be able to map to multiple junctions, especially in cases of short splicing overhangs.

Cheers
Alex

[amolkolte]

saturate the junctions!!

Hi,
Would it be wise to generate the genome with GTF file and re-generate it later with the sample specific SJ.out.tab?

I have more than 100 samples for variant analysis. Another strategy that I thing could be useful is to pull 20-30 samples together (mapping being extremely fast, it was over within 40mins), so that the discovery of novel junctions will get saturated, and then use this final SJ.out.tab to re-generate the genome, and perform sample wise alignment. This has reported ~2,00,000 total junctions, which is quite closer to the number of exons in human genome ~1,80,000

Cheers,
Amol

[dpryan]

Just use the --outFileNamePrefix option to give everything an appropriate prefix.

[dan]

Star is generating files like this in the directory I'm running it from:

• SJ.out.tab
• Log.std.out
• Log.out
• Log.final.out
• Log.progress.out

How can I put those files into a specific directory so that I can run multiple jobs in parallel from a single directory?


Cheers,
Dan.

[serash]

That's great to know! Thanks

[alexdobin]

Hi again,

It makes sense I need to rebuild the genome first and can than use it again for all samples. However, I noticed that rebuilding the genome takes longer than the alignment itself. If I run the 2-pass alignment with --twopass1readsN -1 the outcome is achieved much faster by adding the SJ.out.tab in memory from the first run.
This said, I was wondering if its possible to do the same as the 2-pass run in shared memory, to add the data (give the file by --sjdbFileChrStartEnd?) to an existing genome in memory to avoid rebuilding the entire genome and thus saving time.

cheers,
Dries

This is a good suggestion, and it's already on my (short) TODO-list - adding annotated junctions on the fly, so that there no need to re-generate the whole genome index when you only want to change annotations.

Cheers
Alex

[serash]

Hi again,

It makes sense I need to rebuild the genome first and can than use it again for all samples. However, I noticed that rebuilding the genome takes longer than the alignment itself. If I run the 2-pass alignment with --twopass1readsN -1 the outcome is achieved much faster by adding the SJ.out.tab in memory from the first run.
This said, I was wondering if its possible to do the same as the 2-pass run in shared memory, to add the data (give the file by --sjdbFileChrStartEnd?) to an existing genome in memory to avoid rebuilding the entire genome and thus saving time.

cheers,
Dries

[alexdobin]

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!

Hi,

the --twopass1readsN cannot be used with --genomeLoad LoadAndKeep option, you will need to re-load the genome for every sample. For each of the samples the 2-nd pass genome indices will be different, since the junctions discovered in the 1st pass will change from sample to sample.

For multiple samples, one option is to map (2-pass with --twopass1readsN -1 which will map all reads in the first pass) all the samples together, and use the Read Group tags to mark different samples.

Another option is the manual 2-step operation: run the 1st pass on all samples, collect the detected junctions from all samples, re-generate the genome, and run 2nd pass with the new genome, which can now be loaded with --genomeLoad LoadAndKeep option.

Cheers
Alex