STAR runs slowly

[Babak A]

Hi Alex,

I ran STAR on a couple of thousand experiments and collected around 6.7M junctions. For each junction, I computed the aggregate number of reads in all experiments. Then I filtered them as follows:

1. If annotated, keep it.

2. If unannotated and canonical, keep if has at least 5 aggregate reads in all experiments

3. If unannotated and non-canonical, keep if it has at least 10 aggregate reads.

After filtration and taking the union with GENCODE junctions, 2.3M junctions remained. I generated the index with these junctions.

The problem is that STAR now runs super slowly. On average my experiments have 60M paired end reads, and STAR maps them in 180 seconds. But now, it is 20 times slower:

Compare this:

Jul 19 17:13:15 ..... Started STAR run

Jul 19 17:13:54 ..... Started mapping

Jul 19 18:12:11 ..... Finished successfully

To this:

Jul 08 15:29:09 ..... Started STAR run

Jul 08 15:30:05 ..... Started mapping

Jul 08 15:33:19 ..... Finished successfully

Any suggestions?

Thanks,

Babak

[Alexander Dobin]

Hi Babak,

I would expect some slowdown with a large number of junctions, however, what you observed is very dramatic. Would it be possible for you to share the filtered list of junctions with me, as well as your Log.out file so that I can match your parameters. I will try to figure out the cause of the problem.

Cheers

Alex

[Babak A]

Hi Alex,

Thanks for your response. Attached, please find the log file. The SJDB file is huge, I will find a place to upload it to.

On a side note, I have 128 GB of RAM and a very fast disk ( > 500 MB/s). If I can change any parameters to make STAR run faster, I would appreciate it. Moreover, what would you suggest as the fastest way of aligning, converting to BAM, and then sorting it?

Cheers,

Babak

[Alexander Dobin]

Hi Babak,

I have generated the genome with your 2.2M junctions, and mapped one of my 2x76 datasets.  Mapping speed does indeed drop, although not as dramatically in your runs.

With Gencode 346k junctions, the speed is 537M reads/hour, and with 2.2M junctions it drops to 172M reads/hour, ~3 fold reduction compared to ~20 fold that you observed.

I run it on 12 physical cores, and you are using 64 threads - could you try it with fewer threads, say 16 and 32? It could be that using too many threads leads to degraded performance.

I will look into the speed drop - even a factor of 3 is too much of a drop.

Cheers

Alex

[Babak A]

Hi Alex,

Thanks for taking the time and testing my junctions. 

Please take a look at the two plots I have made using the info in Log.progress.out in two different runs of STAR. I think STAR behaves strangely.

Cheers

Babak

[Babak A]

Sorry, I forgot the X label in the second plot:

[Alexander Dobin]

Hi Babak,

this is indeed weird, the mapping speed should not vary so much in one run. For the first ~65% of the reads, the difference in speed between "2.2M" and "350k" junctions is relatively small - comparable to my observations. Something happens with the last 35% of the reads. The question is whether it's related somehow to the dataset - e.g. much poorer read quality can reduce the mapping rate, or it's some kind of a runtime problem - e.g. memory leakage.

To answer this question, could you please try to map just the last 20% of your reads cutting them out of the .fastq file?

Cheers

Alex

[Babak A]

Hi Alex,

I tried this index on ~80 experiments with various qualities (50%-95% mapped reads) and the average mapping time is 1 hour.

I ran the STAR with the top 15M reads and it did not change much, here is the log:

Time Speed Mapped

14:38:19 634.5 11104066

14:39:37 317.7 12444068

14:40:46 224.9 13118812

14:41:57 174.9 13655622

14:44:09 121.3 13921296

14:45:32 104.0 14325664

14:48:12 80.1 14593648

14:50:36 66.3 14727640

14:55:44 48.7 15000000

Babak

[Babak A]

And when I use only 16 threads instaed of 64 threads, I got the same (slow) running time, but a drastically different pattern:

15:25:23 30.8 538360

15:26:24 39.0 1344704

15:27:29 41.0 2151048

15:28:31 65.5 4567688

15:29:33 109.4 9511104

15:30:54 124.8 13657932

15:32:45 99.3 13928308

15:33:52 89.3 14193982

15:35:50 75.4 14459656

15:38:39 61.7 14727640

15:42:47 48.8 15000000

Babak

[Babak A]

Hi Alex,

I guess I need to run STAR using Gencode junctions, as the larger junction set is not feasible.

Babak

[Alexander Dobin]

Hi Babak,

I am trying to fix the slowdown problem with the large number of junctions. However, I am not sure if it will solve the slowdown at the end of the run problem that you observed - I cannot reproduce that one on my machine unfortunately. Would it be possible for you to share those 15M reads - just in case it's dataset dependent (though it seems unlikely to me).

Cheers

Alex

[Babak A]

Hi Alex,

Thanks for your response.

Unfortunately, I cannot share the data but most likely you already have access to it: SRR595926

I tried that index on several experiments with different percentage of aligned reads and they take one hour to run. I made smaller junction files with 1.5M and 600K junctions and run time nearly linearly scales down.

Cheers,

Babak

[Babak A]

To be more specific:

junctions  running time

350K  ~2m30s 
650K  ~3m30s

1.6M   ~22m

2.2M   ~60m

Babak

[Babak A]

Hi Alex,

I compared the Log.final.out files and the results are quite interesting. As I increase the number of junctions, number of uniquely mapped reads goes down and number of multi-mappers goes up (all in percentile):

Juncs Unique Multi Unmapped

350K 86.81 4.00 9.15

650K 85.66 5.17 9.24

1.1M 83.79 6.91 9.26

1.6M 82.42 8.13 9.41

2.2M 81.06 9.24 9.62

So, now the question is whether adding more junction helps or fires back.

Cheers,

Babak

[Alexander Dobin]

Hi Babak,

I do not have access to GTEX datasets, unfortunately. On my standard testing dataset I can see an increase in % of multimappers at the expense of unique mappers, though not as large as the one you observed:

sjdb     Unique     Multi

350k     89.29%     6.00% 

2.3M     87.49%     7.71%

I looked at the reads that become multi-mappers with 2.3M sjdb, and it appears that most of them map to junctions on the mitochondrial genome.

So I generated genome excluding chrM junctions from your 2.3M set, and got:

sjdb            Unique     Multi 

2.3M, no chrM             88.78%           6.50%

These number are much closer to the 350k sjdb case. Also, the mapping speed in this case is two-fold higher than the 2.3M with chrM case and only ~30% less than the 350k case.

Thus I conclude that at least on my dataset both the mapping slowdown and reduction in the % of unique mappers are caused in large by the sjdb junctions on chrM.

I would suggest that you generate your genome without the junctions on chrM and  try to map your samples to it, hopefully the problems will be resolved as well.

I think it is to be expected that some unique mappers become multi-mappers as you add more and more sjdb junctions, since this effectively adds more possibilities for the reads to align. Note, that by default the --sjdbScore = 2, which means that STAR will try to map aggressively to the sjdb junctions, preferring spliced alignment with 1 mismatch to an unspliced alignment without mismatches. You may want to try to reduce this parameter, though it will lead to yet another slight decrease in the % of unique mappers.

Cheers

Alex

[Babak A]

Hi Alex,

Thank you so very much! As you correctly concluded, chromosome M was the culprit. I removed the junctions in ChrM (<10K) and everything changed dramatically:

1.6M juncs: 7mins (from 22 mins!) - multi-mappers: 4.26% - unique: 86.54%

1.1M juncs: 4.5mins - multi-mappers:  4.17% - unique: 86.64%

I cannot believe what less than 10K junctions can do.

Cheers,

Babak