Fastq file from SRA database not mapping to to the genome

Question

[amoltej]

Hi I have downloaded SRA file from the database and generated fastq files using SRA-toolkit. I tried to map these to custom genome file (fasta seq) using STAR, Bowtie, Bowtie2, BWA. But none of them worked.

I tried to download fastq files from EBI website but it didn't work. When I opened the fastq file I see weird encoding when I compared it to the fastq file generated by my own experiment. please see below - SRA fastq file-

@SRR867425.1 TEST690_0001:2:1:1013:1620/1
NACCACTATTGGACAAATCCAGGGAACATNGTCACTTCAGAACCAGAGTGCTTTTATTAGGCTTTTAGTCTGCTGTCCAGCTCTC
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%@
@SRR867425.2 TEST690_0001:2:1:1015:20229/1
NTGCAATTCCGGCGAAGGAACCACTGTTTCAGGTCGAGTAAGCACGAATCCTTCTTCCAGTAACGAGACCCCTGAAAATCCTCCC
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%A
@SRR867425.3 TEST690_0001:2:1:1015:8397/1
NAAACACGTAATATGATGCTTCTACATAGCCATCTAACTGAGTCAGTTGAGACCATAAACCACTGCTTCTGCCATGAGGCAGGAA
+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%=

Fastq file from my own experiment

@HWI-700666F:84:C7GUUANXX:4:1101:2067:1954 1:N:0:AGTCAA
ATTTTTATTTTAAAAAGAGGATGTGTTTCCAATCAGTCTTTGCGCTTCTTCGATCCTTGTTCCTTTTGTGCTAAAATAATTAGCAGCCTTCTGAAGATG
+
=30:FFBFFGGG1CD@GGE>FGG@==:FGGGGGBGCGGGGG1@CGGGGGGG0CDGGGCDEGE>FGGGGGGFEG11<DCF>G:11?:FB0BDF:000BCB
@HWI-700666F:84:C7GUUANXX:4:1101:2004:1973 1:N:0:AGTCAA
TGGCCTTTATGAGTGGTCAAAGCTGCCTTGAGGAGAATTGTTTAATTCATCAACCTCCTGTATAACAGATTCATCTTCCAAACGCTGTCGATCAATTAA
+
<<ABFGGGGGGGGGGGGGGGGGGGDCFGG>DGGGDGGGCGGGGG1FGEGGGGGGGGGGFGGGGEGGGG1=EFGGGG1FGGCFGBEF=FGBGGGGGGEED

Can someone please help me to sort the issue? Thanks in advance

[Nitin Narwade]

What was the error? Can please post it.

[amoltej]

When I ran STAR this is the last part of the log.out file -

Writing 7077888 bytes into ./Genome ; empty space on disk = 356075089920 bytes ... done
SA size in bytes: 238552
Jul 11 16:20:35 ... writing Suffix Array to disk ...
Writing 238552 bytes into ./SA ; empty space on disk = 356068012032 bytes ... done
Jul 11 16:20:35 ... writing SAindex to disk
Writing 8 bytes into ./SAindex ; empty space on disk = 356067770368 bytes ... done
Writing 120 bytes into ./SAindex ; empty space on disk = 356067770368 bytes ... done
Writing 1565873491 bytes into ./SAindex ; empty space on disk = 356067770368 bytes ... done
Jul 11 16:20:37 ..... Finished successfully
DONE: Genome generation, EXITING`

[GenoMax]

That log appears to be just for generation of genome indexes. Have you done the actual alignment? Please post the full STAR command lines you have used for genome index generation and for alignment runs.

[amoltej]

first I ran this STAR --runMode genomeGenerate --genomeDir . --genomeFastaFiles MYB_and_MYC.fasta --runThreadN 5

Then STAR --runThreadN 14 --sjdbOverhang 50 --readFilesCommand zcat --sjdbGTFfile myb_myc.gtf --readFilesIn SRR867425_1.fastq.gz SRR867425_2.fastq.gz --runMode alignReads --outReadsUnmapped Fastx --limitBAMsortRAM 8729257684 --genomeDir . --outSAMtype BAM SortedByCoordinate --outSAMstrandField intronMotif --outFilterIntronMotifs RemoveNoncanonical --twopassMode Basic --outSAMattrRGline ID:"MYB_MYC" --outFileNamePrefix "MYB_MYC"_ --quantMode GeneCounts

[GenoMax]

I suggest that you create a separate directory to hold the STAR indexes and supply it to --genomeDir option (instead of just ./).

[amoltej]

its still not working. exiting with the error message Segmentation fault (core dumped) at the 1st pass mapping

[GenoMax]

At least we have made some progress. How much memory are you assigning to this job (or should ask how much memory you have available)?

[amoltej]

I am not assigning any particular number but I am running with 15 threads with 64GB ram. No other program is running on machine

[GenoMax]

Try reducing the number to 4 and see if that prevents the seg fault.

Answer 1

[GenoMax]

% is a valid value in Sanger fastq format. It is not unusual to see some reads at the beginning of the files (at edge of flowcells) have low Q scores like this. I don't think this is causing the alignment issue.

map these to custom genome file (fasta seq)

Can you clarify what this is? A subset fasta formatted file you used as a custom genome?

Edit: Three example reads posted in the example above are matching Citrus species by blast. That matches what the SRR867425 submission.

Edit 2: Reads from your experiment are from an aphid so you must be looking at a host-pest experiment of some sort. Are you trying to separate reads from these species?

[amoltej]

Thanks for reply. Yes, both files are from different organisms. I just used them to compare if there is any issue with the fastq file format. Here is the fasta file part that i am using. this is basically few genes I selected for my study -

>NM_001112236_1_Zea_mays_colored_plant_1_b1
ATGGCGCTCTCAGCTTCCCCGGCTCAGGAAGAACTGCTGCAGCCTGCTGGGAGGCCGTTGAGGAAGCAGCTT
GCTGCAGCCGCGAGGAGCATCAACTGGAGCTATGCCCTCTTCTGGTCCATTTCAAGCACTCAACGACCTC
GGGTGCTGACGTGGACGGACGGGTTCTACAATGGCGAGGTGAAGACGCGTAAGATCTCCCACTCCGTGGA
GCTGACAGCCGACCAGCTGCTCATGCAGAGGAGCGAGCAGCTCCGGGAGCTCTACGAGGCCCTCCGGTCC
GGCGAGTGCGACCGCCGCGGCGCGCGGCCGGTGGGCTCGCTGTCGCCGGAGGACCTCGGGGACACCGAGT
>DQ275529_1_Antirrhinum_majus_ROSEA1
ATGGAAAAGAATTGTCGTGGAGTGAGAAAAGGTACTTGGACCAAAGAAGAAGACACTCTCTTGAGGCAAT
GTATAGAAGAGTATGGTGAAGGGAAATGGCATCAAGTTCCACACAGAGCAGGGTTGAACCGGTGTAGGAA
GAGTTGCAGGCTGAGGTGGTTGAATTATCTGAGGCCAAATATCAAAAGAGGTCGGTTTTCGAGAGATGAA
GTGGACCTAATTGTGAGGCTTCATAAGCTGTTGGGTAACAAATGGTCGCTGATTGCTGGTAGAATTCCTG
GAAGGACAGCTAATGACGTGAAGAACTTTTGGAATACTCATGTGGG

[GenoMax]

Nothing you are doing seems to be out of ordinary (assuming you created the indexes from the fasta file properly etc). Perhaps your data does not have reads that map to these regions. Have you considered that possibility?

[amoltej]

the index is creating successfully. I considered the possibility of the no mapped reads. but it should at least give me a report or something with the STAR program. or empty BAM file. There is nothing like that. And that is why its confusing to me!

[GenoMax]

1. I made a fake read including sequence from your reference (third read).

   @HWI-700666F:84:C7GUUANXX:4:1101:2067:1954 1:N:0:AGTCAA ATTTTTATTTTAAAAAGAGGATGTGTTTCCAATCAGTCTTTGCGCTTCTTCGATCCTTGTTCCTTTTGTGCTAAAATAATTAGCAGCCTTCTGAAGATG + =30:FFBFFGGG1CD@GGE>FGG@==:FGGGGGBGCGGGGG1@CGGGGGGG0CDGGGCDEGE>FGGGGGGFEG11<dcf>G:11?:FB0BDF:000BCB @HWI-700666F:84:C7GUUANXX:4:1101:2004:1973 1:N:0:AGTCAA TGGCCTTTATGAGTGGTCAAAGCTGCCTTGAGGAGAATTGTTTAATTCATCAACCTCCTGTATAACAGATTCATCTTCCAAACGCTGTCGATCAATTAA + <<abfgggggggggggggggggggdcfgg>DGGGDGGGCGGGGG1FGEGGGGGGGGGGFGGGGEGGGG1=EFGGGG1FGGCFGBEF=FGBGGGGGGEED @HWI-700666F:84:C7GUUANXX:4:1101:2005:1974 1:N:0:AGTCAA GGTGAAGGGAAATGGCATCAAGTTCCACACAGAGCAGGGTTGAACCGGTGTAGGAAGAGTTGCAGGCTGAGGTGGTTGAATTATCTGAGGCCAAATATC + <<abfgggggggggggggggggggdcfgg>DGGGDGGGCGGGGG1FGEGGGGGGGGGGFGGGGEGGGG1=EFGGGG1FGGCFGBEF=FGBGGGGGGEED

2. Put your reference in gen.fa (fasta sequence above).

3. Aligned the fastq file to reference using bbmap.sh. I can see the third spiked read aligning properly in the resulting sam file.

   HWI-700666F:84:C7GUUANXX:4:1101:2067:1954 1:N:0:AGTCAA 4 * 0 0 * * 0 0 ATTTTTATTTTAAAAAGAGGATGTGTTTCCAATCAGTCTTTGCGCTTCTTCGATCCTTGTTCCTTTTGTGCTAAAATAATTAGCAGCCTTCTGAAGATG =30:FFBFFGGG1CD@GGE>FGG@==:FGGGGGBGCGGGGG1@CGGGGGGG0CDGGGCDEGE>FGGGGGGFEG11<dcf>G:11?:FB0BDF:000BCB HWI-700666F:84:C7GUUANXX:4:1101:2004:1973 1:N:0:AGTCAA 4 * 0 0 * * 0 0 TGGCCTTTATGAGTGGTCAAAGCTGCCTTGAGGAGAATTGTTTAATTCATCAACCTCCTGTATAACAGATTCATCTTCCAAACGCTGTCGATCAATTAA <<abfgggggggggggggggggggdcfgg>DGGGDGGGCGGGGG1FGEGGGGGGGGGGFGGGGEGGGG1=EFGGGG1FGGCFGBEF=FGBGGGGGGEED HWI-700666F:84:C7GUUANXX:4:1101:2005:1974 1:N:0:AGTCAA 0 DQ275529_1_Antirrhinum_majus_ROSEA1 85 44 99= * 0 GGTGAAGGGAAATGGCATCAAGTTCCACACAGAGCAGGGTTGAACCGGTGTAGGAAGAGTTGCAGGCTGAGGTGGTTGAATTATCTGAGGCCAAATATC <<abfgggggggggggggggggggdcfgg>DGGGDGGGCGGGGG1FGEGGGGGGGGGGFGGGGEGGGG1=EFGGGG1FGGCFGBEF=FGBGGGGGGEED NM:i:0 AM:i:44

[amoltej]

I was wondering if you use the fastq file with % sign for making fake read and then mapping would it work?

[GenoMax]

Those reads with % quality scores are poor quality and will fail to map. Not all reads in the SRA file have that type of Q scores, correct?

[amoltej]

It is there in almost all the reads. but the number of % signs are different in different reads.