RNA-STAR, mapping problem (number of reads)
0
0
Entering edit mode
7.0 years ago
mohikran • 0

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 • 7.1k views ADD COMMENT 1 Entering edit mode You should not put a space between the -- and the parameter name. ADD REPLY 0 Entering edit mode What is written in the Log.final.out ? ADD REPLY 0 Entering edit mode 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%  ADD REPLY 2 Entering edit mode Try that for the mapping: STAR --runThreadN 5 --genomeDir$pathINDEX --readFilesIn ./MYFILE.fastq

0
Entering edit mode

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

0
Entering edit mode

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
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
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%
% of reads unmapped: too many mismatches |       0.00%
% of reads unmapped: too short |       20.44%
% of reads unmapped: other |       1.44%
Number of chimeric reads |       0
% of chimeric reads |       0.00%

0
Entering edit mode

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