This is a test version of Biostars. For the public version, visit https://www.biostars.org.
RNA-STAR, mapping problem (number of reads)

Hello everyone,

I have some RNA-seq (Illumina/ single-end) data to analyze, provided by my boss and I have one problem with the mapping using RNA-STAR.

So far, I am able to use Fastqc, some trimming, generating my genome, but I have a problem with the mapping step. I would like to get the number of read for each gene (and then I want to use Deseq2 with R), but I don't have any file created with these information after running my script, and my Log.out file does not really help me this time. I have the following files created: Log.final.out, Log.out and Log.progress but not any BAM/SAM file.

When I have a look in the log, it is said that the mapping is done, but nothing else...

I have to admit that I am a bit lost and if someone could help me that would be great. I guess that the problem comes from the generation of the genome (a wrong annotation file?) or one of my parameter during the mapping.

Thank you and have a nice day! :)

This is how I generated my *genome *(All the files come from ensembl, genome GRCm38):

#===== BUILD GENOME INDEX
cd $pathINDEX
cp MOUSE_index_GRCm38.fa $SNIC_TMP #This one contains the merged fasta sequences of the chromosomes.
cd $pathGTF
cp Mus_musculus.GRCm38.83.chr.gtf.gz $SNIC_TMP #The file containing the annotations
cd $SNIC_TMP
gzip -d Mus_musculus.GRCm38.83.chr.gtf.gz
STAR --runMode genomeGenerate --genomeDir $pathFINAL --runThreadN 5 --sjdbGTFfile Mus_musculus.GRCm38.83.ch
echo "########## END ##########"

And the mapping:

#===MAPPING
gzip -d Sca1_CD51_ko_131230_CACCTC_L003_R1_001.fastq.gz
STAR -- runThreadN 5 -- genomeDir $pathINDEX -- readFilesIn ./ MYFILE.fastq
echo "########## END ##########"

What is written in Log.out :

Finished loading and checking parameters
Reading genome generation parameters:
versionGenome                 20201        ~RE-DEFINED
genomeFastaFiles              MOUSE_index_GRCm38.fa        ~RE-DEFINED
genomeSAindexNbases           14     ~RE-DEFINED
genomeChrBinNbits             18     ~RE-DEFINED
genomeSAsparseD               1     ~RE-DEFINED
sjdbOverhang                  100     ~RE-DEFINED
sjdbFileChrStartEnd           -        ~RE-DEFINED
sjdbGTFfile                   Mus_musculus.GRCm38.83.chr.gtf     ~RE-DEFINED
sjdbGTFchrPrefix              -     ~RE-DEFINED
sjdbGTFfeatureExon            exon     ~RE-DEFINED
sjdbGTFtagExonParentTranscripttranscript_id     ~RE-DEFINED
sjdbGTFtagExonParentGene      gene_id     ~RE-DEFINED
sjdbInsertSave                Basic     ~RE-DEFINED
Genome version is compatible with current STAR version
Number of real (reference) chromosmes= 21
1       10      130694993       0
2       11      122082543       130809856
3       12      120129022       252968960
4       13      120421639       373293056
5       14      124902244       493879296
6       15      104043685       618921984
7       16      98207768        722993152
8       17      94987271        821297152
9       18      90702639        916455424
10      19      61431566        1007419392
11      1       195471971       1069023232
12      2       182113224       1264582656
13      3       160039680       1446772736
14      4       156508116       1606942720
15      5       151834684       1763704832
16      6       149736546       1915748352
17      7       145441459       2065694720
18      8       129401213       2211184640
19      9       124595110       2340683776
20      X       171031299       2465464320
21      Y       91744705        2636644352
--sjdbOverhang = 100 taken from the generated genome
Started loading the genome: Fri Feb 19 18:37:29 2016

checking Genome sizefile size: 2779951453 bytes; state: good=1 eof=0 fail=0 bad=0
checking SA sizefile size: 22265276075 bytes; state: good=1 eof=0 fail=0 bad=0
checking /SAindex sizefile size: 1565873619 bytes; state: good=1 eof=0 fail=0 bad=0
Read from SAindex: genomeSAindexNbases=14  nSAi=357913940
nGenome=2779951453;  nSAbyte=22265276075
GstrandBit=32   SA number of indices=5397642684
Shared memory is not used for genomes. Allocated a private copy of the genome.
Genome file size: 2779951453 bytes; state: good=1 eof=0 fail=0 bad=0
Loading Genome ... done! state: good=1 eof=0 fail=0 bad=0; loaded 2779951453 bytes
SA file size: 22265276075 bytes; state: good=1 eof=0 fail=0 bad=0
Loading SA ... done! state: good=1 eof=0 fail=0 bad=0; loaded 22265276075 bytes
rna-seq
• 8,842 views • link
updated by Ram Ram • written by mohikran • mohikran

Hello, and first thank you for your answer :) .

I will do as you recommend, removing the spaces after the "--" and trying to use featureCounts.

This is what I've got in my Log.final.out. It seems that the mapping does not work, that's why I am wondering if my genome is properly generated (wrong GTF file?).

                                 Started job on |       Feb 19 18:37:28
                             Started mapping on |       Feb 19 18:43:38
                                    Finished on |       Feb 19 18:43:40
       Mapping speed, Million of reads per hour |       0.00

                          Number of input reads |       0
                      Average input read length |       0
                                    UNIQUE READS:
                   Uniquely mapped reads number |       0
                        Uniquely mapped reads % |       0.00%
                          Average mapped length |       0.00
                       Number of splices: Total |       0
            Number of splices: Annotated (sjdb) |       0
                       Number of splices: GT/AG |       0
                       Number of splices: GC/AG |       0
                       Number of splices: AT/AC |       0
               Number of splices: Non-canonical |       0
                      Mismatch rate per base, % |       -nan%
                         Deletion rate per base |       0.00%
                        Deletion average length |       0.00
                        Insertion rate per base |       0.00%
                       Insertion average length |       0.00
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       0
             % of reads mapped to multiple loci |       0.00%
        Number of reads mapped to too many loci |       0
             % of reads mapped to too many loci |       0.00%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |       0.00%
                 % of reads unmapped: too short |       0.00%
                     % of reads unmapped: other |       0.00%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
                            % of chimeric reads |       0.00%
• 1 views • link
updated by Ram Ram • written by mohikran • mohikran

Try that for the mapping:

STAR --runThreadN 5 --genomeDir $pathINDEX --readFilesIn ./MYFILE.fastq
• 1 views • link
updated by Ram Ram • written by Nicolas Rosewick Nicolas Rosewick

And to have a count table you should use featureCounts on your sam/bam file generated by STAR

• 0 views • link
updated by Ram Ram • written by Nicolas Rosewick Nicolas Rosewick

Hello,

It workded, now I have a out.sam file created! Thanks again, Now I will use featureCounts as you said on my sam file.

However if I can ask a last thing. I just realized that the pipeline that I am using in not performing any assembly of the reads. My reads are shorts (~50) and directly mapped to the genome. Is it a problem? For my analyses, we're just planning to get the number of reads and use Deseq.

Best regards,

Here is my log file:

                                 Started job on |       Feb 23 13:58:12
                             Started mapping on |       Feb 23 14:03:52
                                    Finished on |       Feb 23 14:09:56
       Mapping speed, Million of reads per hour |       319.89

                          Number of input reads |       32344665
                      Average input read length |       50
                                    UNIQUE READS:
                   Uniquely mapped reads number |       20927395
                        Uniquely mapped reads % |       64.70%
                          Average mapped length |       48.22
                       Number of splices: Total |       981860
            Number of splices: Annotated (sjdb) |       907516
                       Number of splices: GT/AG |       953422
                       Number of splices: GC/AG |       6598
                       Number of splices: AT/AC |       1208
               Number of splices: Non-canonical |       20632
                      Mismatch rate per base, % |       1.13%
                         Deletion rate per base |       0.01%
                        Deletion average length |       1.49
                        Insertion rate per base |       0.01%
                       Insertion average length |       1.24
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       4022138
             % of reads mapped to multiple loci |       12.44%
        Number of reads mapped to too many loci |       317960
             % of reads mapped to too many loci |       0.98%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |       0.00%
                 % of reads unmapped: too short |       20.44%
                     % of reads unmapped: other |       1.44%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
                            % of chimeric reads |       0.00%
• 0 views • link
updated by Ram Ram • written by mohikran • mohikran

With RNA-seq you are not looking to assemble reads into contigs as you would if you were doing genomic sequencing.

• 0 views • link
written by doliv071 • doliv071

You should not put a space between the -- and the parameter name.

• 0 views • link
updated by Ram Ram • written by Nicolas Rosewick Nicolas Rosewick

0 answers

No answers yet.

Log in to answer this question.

More posts like this