[SamCurt]

Just a quick question here. Is the parameters file used with --parametersFile just a list of command-line options in the same way I type in the console?

[Juegos 2 friv 4school]

Thank your article. very helpful article. thank you very much.

[harrike]

Hi Alex,

Thanks for your suggestions.

I manually checked a couple of reads as Genomax suggested, and find the major reason of this low mapping rate is because that most of the reads have adapter, due to the poor construction of RNA-seq library. What I am trying to do is to trim the adapter and do the mapping again. The read quality is good per FastQC check.

I will try to relax --alignEndsType option, and see if the mapping will become better or not.

Rui

[alexdobin]

Hi Rui,

here are a few suggestions in addition to @GenoMax's suggestion.

1. You are using --alignEndsType EndToEnd, which requires end-to-end alignment for each read (no soft clipping). This might be too harsh for longer reads, which are more likely to have poor quality tails, adapters at the ends etc. Please try to map without this option.
2. Map read1 and read2 separately - you may have a problem with one of the reads.
3. Check sequencing quality by plotting quality scores vs position in read (Illumina pipelines typically produce these plots). If sequencing quality drops towards the ends of the reads for a substantial portion of the reads, this would explain poor mappability.

Cheers
Alex

[GenoMax]

I suggest that you start by looking at a few (10-20) unmapped reads and blast them against nt to see what they are aligning to. You may be surprised by what you find and it may provide an explanation for the low % alignment.

[harrike]

Hi Alex,

This time I am using STAR to another set of data, which are strand-specific, paired-end, and of 150 bp read length. The command I used is "

STAR-STAR_2.4.2a/bin/Linux_x86_64/STAR --genomeDir Zmay_AGPv2_STAR_index/ --runThreadN 24 --readFilesIn Zm_ant_02_07a_TGACCA_L001_R1_001.fastq Zm_ant_02_07a_TGACCA_L001_R2_001.fx/ --runThreadN 24 --readFilesIn Zm_TGACCA_L001_R1_001.fastq Zm_TGACCA_L001_R2_001.fastq --outSAype EndToEnd --outFilterIntronMotifs RemoveNoncanonical --outFilterType BySJout --outFileNamePrefix Zm_antMtype BAM Unsorted --outFilterMultimapNmax 20 --alignIntronMax 10000 --alignMatesGapMax 10000 --alignEndsType EndToEnd --outFilterIntronMotifs RemoveNoncanonical --outFilterType BySJout --outFileNamePrefix Zm_TGACCA_L001_R1R2_

There are 56.34 of reads unmapped (short), see the Log.final.out file below:

Started job on | Feb 19 07:56:41
Started mapping on | Feb 19 07:57:00
Finished on | Feb 19 08:02:30
Mapping speed, Million of reads per hour | 209.57

Number of input reads | 19210657
Average input read length | 302
UNIQUE READS:
Uniquely mapped reads number | 5863394
Uniquely mapped reads % | 30.52%
Average mapped length | 301.78
Number of splices: Total | 5527054
Number of splices: Annotated (sjdb) | 5329770
Number of splices: GT/AG | 5445210
Number of splices: GC/AG | 76305
Number of splices: AT/AC | 5539
Number of splices: Non-canonical | 0
Mismatch rate per base, % | 0.83%
Deletion rate per base | 0.07%
Deletion average length | 3.02
Insertion rate per base | 0.07%
Insertion average length | 1.98
MULTI-MAPPING READS:
Number of reads mapped to multiple loci | 1335209
% of reads mapped to multiple loci | 6.95%
Number of reads mapped to too many loci | 5625
% of reads mapped to too many loci | 0.03%
UNMAPPED READS:
% of reads unmapped: too many mismatches | 4.69%
% of reads unmapped: too short | 56.34%
% of reads unmapped: other | 1.46%

What are the possible reason of this low-mapping rate? Thanks,

Rui

[harrike]

Hi Genomax,

Thanks for providing the info. It is quite helpful. I am clear now.

Rui

[GenoMax]

Does the 2nd column indicate the strand? I check a few lines of my data. The frequency and number (2nd column value) are listed below.

Frequency Number
847 147
75 163
902 339
94 355
97 403
909 419
71 83
847 99

419/339 pair stands for the "+" strand, 147/99 the "-" strand"? What others?

See the added info in my post above.

[harrike]

Does the 2nd column indicate the strand? I check a few lines of my data. The frequency and number (2nd column value) are listed below.

Frequency Number
847 147
75 163
902 339
94 355
97 403
909 419
71 83
847 99

419/339 pair stands for the "+" strand, 147/99 the "-" strand"? What others?

[GenoMax]

I am also interested to know what the second column stands for in the above example? I am not able to figure it out.

Thanks,

Rui

Second column is bitwise FLAG. Check section 1.4 from full SAM spec: https://samtools.github.io/hts-specs/SAMv1.pdf

Bitwise FLAG also tells you the strand information (https://www.biostars.org/p/59388/). You can "translate" the FLAG using this tool: https://broadinstitute.github.io/pic...ain-flags.html

[harrike]

Hi Alex,

I have a simple question that is tell the strand of a read mapped to (strand-specific data). I viewed my output .bam file using "sambamba view". Is it possible to tell tell the strand from the example given below?

HWI-ST1395:97:d29b4acxx:5:2105:5630:92677 99 chr4 15003609 3 100M = 15003773 264 GCATAATTATGAAGATATTTTTTTATATGGTCCGACTATATAAAATAAATATCATAACGGAAATCAACTTGAACAATGCTGTACCAAAACTACAGT
GGTT CCCFFFFFHHHHGJIJJIJJJIJJJJJJJJGIJJIIJJJIJJJIJJGGJJIJJJJIJJIJIJHGGGFFDFDFECEEECDDDDDDCDDDDDDCCDD@@C@B NH:i:2 HI:i:1 AS:i:198 nM:i:0

HWI-ST1395:97:d29b4acxx:7:2301:9293:56479 99 chr4 15003609 3 100M = 15003701 192 GCATAATTATGAAGATATTTTTTTATATGGTCCGACTATATAAAATAAATATCATAACGGAAATCAACTTGAACAATGCTGTACCAAAACTACAGT
GGTT CCCFFFFFHHHHHJJJIIJJJJJJHJJJJJJJIIJJJJIIJJGJIIJJJJIJJJJGGIJIGJFEHHEFFFEFECEEECDDDACCCDCDDBDDCDD>@CCD NH:i:2 HI:i:1 AS:i:198 nM:i:0

If not, how can I do that (get the strand information). I am also interested to know what the second column stands for in the above example? I am not able to figure it out.

Thanks,

Rui

[alexdobin]

I am attempting to crossover to STAR from TopHat for my RNAseq workflow but I had a question regarding its capabilities.

Specifically, TopHat had an option to include mate-inner-dist and mate-std-dev but relied on RSeQC to input these values for each unique samples.

I was wondering if STAR provides an option to streamline/automate this process. Thank you very much for the awesome tool!

Hi @umn_bist

STAR does not need or use the "expected" insert length information for mapping - so yo do not have to worry about it.
For mapping of DNA-seq data we can decide whether an alignment has insert size that's too big or too small by comparing it "expected" distribution of insert size.
However, in RNA-seq the unsequenced portion of insert between the mates can contain a splice junctions.
This means that we cannot simply calculate the insert size from an alignment to compare it with "expected" insert size.

Cheers
Alex

[emixaM]

Yes I am using a schedule (MOAB) but never had any problem accessing all the RAM of the nodes.

[GenoMax]

Are you using a job scheduler then? Perhaps the scheduler requires an explicit request for more RAM when not using default (e.g. with LSF it would be -M 20).

[emixaM]

Yes, I am the only one, it runs on a compute node of a HPC.

Here is the top of a compute node (no STAR running):

$ top -M
top - 14:33:48 up 3 days, 23:33,  0 users,  load average: 0.00, 1.98, 5.10
Tasks: 333 total,   1 running, 332 sleeping,   0 stopped,   0 zombie
Cpu(s):  0.0%us,  0.0%sy,  0.0%ni, 99.9%id,  0.0%wa,  0.0%hi,  0.0%si,  0.0%st
Mem:    23.584G total, 1948.043M used,   21.682G free,    0.000k buffers
Swap:    0.000k total,    0.000k used,    0.000k free, 1073.691M cached

  PID USER      PR  NI  VIRT  RES  SHR S %CPU %MEM    TIME+  COMMAND
31824 emixaM   20   0 27620 1620 1072 R  0.3  0.0   0:00.06 top
...