STAR mapping pipeline with 2-pass for multiple samples

[Rajesh Detroja]

Dear Alex,

I am running an STAR 2.7.0f mapping pipeline with 2-pass mode for multiple samples of different diseases and healthy samples as follows:

1) Indexing genome with annotations

STAR --runMode genomeGenerate --genomeDir ~/db/hg38/ --genomeFastaFiles ~/db/hg38/hg38.fa --sjdbGTFfile ~/db/hg38/hg38.gtf --runThreadN 30 --sjdbOverhang 89

Note:

- Indexing for maximum read lengh 90 bp.

2) 1-pass mapping with indexed genome

STAR --genomeDir ~/db/hg38/ --readFilesIn sample1.R1.fastq.gz sample1.R2.fastq.gz --readFilesCommand zcat --outSAMunmapped Within --outFileNamePrefix sample1. --runThreadN 30

Notes:

- The same command has been run for multiple samples in the for loop, therefore, it will generate SJ.out.tab file for each sample.

- Next, I have copied SJ.out.tab files of all the samples into the a single folder "SJ_out"

3) Indexing genome with annotations and SJ.out.tab files

STAR --runMode genomeGenerate --genomeDir ~/db/hg38/SJ_Index/ --genomeFastaFiles ~/db/hg38/SJ_Index/hg38.fa --sjdbGTFfile ~/db/hg38/SJ_Index/hg38.gtf --runThreadN 30 --sjdbOverhang 89 --sjdbFileChrStartEnd SJ_out/*.SJ.out.tab

Note:

- Again indexing for maximum read lengh 90 bp.

2) 2-pass mapping with new indexed genome with annotations and SJ.out.tab files

STAR --genomeDir ~/db/hg38/SJ_Index/ --readFilesIn sample1.R1.fastq.gz sample1.R2.fastq.gz --readFilesCommand zcat --outSAMunmapped Within --outFileNamePrefix sample1. --runThreadN 30

Notes:

- Again, same command has been run for multiple samples in the for loop, therefore, it will genreated mapping files for each sample.


Could you please help me to check if all the commands I am running correctly or do you have any suggestions?

[Alexander Dobin]

Hi Rajesh,

the commands look good to me. The only concern is that, if you have too many samples, and too many novel junctions are detected across all samples, the 2nd pass may become too slow and the number of multimappers may be increased.

Then you would need to filter the set of junctions you are using for the 2-pass genome.

Cheers

Alex

[Rajesh Detroja]

Dear Alex,

Thanks for the confirming the commands.

According to one of your recent post you have suggested following possible criterias to filter the merged junction file of multiple samples, Is that correct?

1. Filter out the junctions on chrM, those are most likely to be false. That's mean remove all chrM records.

2. Filter out non-canonical junctions (column5 == 0). That's mean (column5 > 0).

3. Filter out junctions supported by multi mappers only (column7==0). That's mean (column7 > 0)

4. Filter out junctions supported by too few reads (e.g. column7<=2). That's mean (column7 > 2)


One more question, I should merge junction files of diseases and healthy samples into one or I need to make separate index for diseases and healhy for 2-pass?

[Alexander Dobin]

Hi Rajesh,

this looks like a good filtering strategy. However, you may need to make filtering more stringent if you see a significant slowdown of the 2nd pass, or increase in multimappers.

If you are planning to compare novel splice junctions between healthy and diseased samples, it's better to combine the junctions from both healthy and diseased samples into one index.

Cheers

Alex

[Rajesh Detroja]

Hello Alex,

Thanks for your help.

I have calculated the number of junctions after filtering as follows and I have decided to take filtration with 340,968 Junctions remained:

What do yout think?

cat *.tab | awk '($5 > 0 && $7 > 2)' | wc -l : 32,985,210 Junctions remained


cat *.tab | awk '($5 > 0 && $7 > 5)' | wc -l: 26,981,895 Junctions remained


cat *.tab | awk '($5 > 5 && $7 > 5)' | wc -l: 15,359 Junctions remained

cat *.tab | awk '($5 > 2 && $7 > 5)' | wc -l: 243,547 Junctions remained

cat *.tab | awk '($5 > 1 && $7 > 5)' | wc -l: 13,645,142 Junctions remained

cat *.tab | awk '($5 > 1 && $7 > 5)' | grep -v "MT" | wc -l: 13,645,060 Junctions remained

cat *.tab | awk '($5 > 2 && $7 > 2)' | wc -l: 340,968 Junctions remained

[Alexander Dobin]

Hi Rajesh,

I think you are counting multiple times the same junctions from many samples. You would need to "collapse" the junctions, i.e. remove the duplicate ones after filtering, e.g.

cat *.tab | awk '($5 > 0 && $7 > 2 && $6==0)' | cut -f1-6 | sort | uniq |wc -l

I added $6==0 condition to avoid counting annotated junctions since they are always included from the GTF.

You can feed the resulting small file to STAR, or the original non-collapsed file - STAR will collapse the junctions before inserting them.

Cheers

Alex

[Rajesh Detroja]

Dear Alex,

Thank you so much for your help!

Finally, I have total 501,404 junctions remaining that I will feed to START for generating a new index for 2-pass mapping.

cat *.tab | awk '($5 > 0 && $7 > 2 && $6==0)' | cut -f1-6 | sort | uniq | wc -l

I got total 501,404 records.

Moreover, 2-pass mapping is not slow at all, I feel as much fast as 1-pass.

Please mark this post as completed, I hope it this post will be useful for many peoples too.

[Alexander Dobin]

Hi Rajesh!

Thanks a lot for sharing your experience!

It will definitely be helpful to other users.

Cheers

Alex

[Pedro Milanez-Almeida]

Dear Dr. Dobin and Dr. Detroja,

If I want to keep the resulting file after using 

cat *.tab | awk '($5 > 0 && $7 > 2 && $6==0)' | cut -f1-6 | sort | uniq |wc -l

all I need to do is replace "| wc -l" in the end by "> SJ.filtered", right?

Thank you very much!

Pedro

PS: using STAR2.4.2a

[Rajesh Detroja]

Dear Pedro,

Yes, you are correct! You can keep resulting file as follows:

cat *.tab | awk '($5 > 0 && $7 > 2 && $6==0)' | cut -f1-6 | sort | uniq > SJ.filtered.tab

[Pedro Milanez-Almeida]

Awesome! Thank you very much! Also thanks for sharing your code and explaining what it means. It's very helpful.

[Sreya Mukherjee]

Hello

I am quite new to STAR hence was wondering if you could share the loop part of your code to take care of multiple samples being mapped in STAR ?

[rohit satyam]

Hi, I wanted to know if removing the mitochondrial junctions should be practiced for all genomes such as the rice genome too?

[Alexander Dobin]

Hi Rohit,

yes, it makes sense to remove mitochondrial splice junctions for all species.

Cheers

Alex

[Bekah W]

Hi Alex,

Is there a reason why non-canonical splice sites should be removed for running the second pass?

Best wishes,
Bekah

[Alexander Dobin]

Hi Bekah,

in general, non-canonical junctions have a much higher rate of false positives, and so their detection is unreliable.

If you are not specifically interesting in non-canonical junctions, it's better to filter them out before the 2nd pass.

Cheers

Alex