some help ! BAM to Fastq

[elitatoun]

Dear all

I'm new on your forum and on the NGS bio informatic analysis !

Sorry for my english !

I need some help from specialists...i try to analyse NGS results for small non coding RNA...

until now i have made :
checked the qulity of the reads
removed the adaptators...
launched my first map (bowtie). And i have found that the half of my reads don't mapp my genom of reference...

and i try to understand why !

so i try to keep all the unmapped reads, to see what they are :
i have used bwa to keep my unmapped reads and i have obtained a SAI file.
I have converted the SAI file into BAM file

Now i will try to creat a new BAM file with only unmapped reads with this Command :

samtools view fichier bam.bam | awk '$3 =="*"' | samtools view -bS - > no_maps.bam

and after i will convert this BAM file into Fastq file with :


java -jar -Xmx6g ~/bin/SamToFastq.jar INPUT=no_maps.bam F=no_map.fastq VALIDATION_STRINGENCY=LENIENT

A very important think :
i have found these 2 commands, but i don't understand everything...that why i need your help !

If someone can explain me, it will be great

Thx

[GenoMax]

I am moving this post to a new thread for better visibility.

@elitatoun: If you were only interested in capturing reads that did not map with bowtie you could have done that using the "--un filename" option. This would have captured all reads that did not align to a new file.

Are you trying to find what these reads are? You could take a few of them (you will need to convert them to fasta format) and blast them at NCBI to get a general idea of what they are.

[elitatoun]

@GenoMax
thank you for your reply ! i will try also your solution thx
and yes i'm trying to find what these reads are....
yes i will blast (a small part of them) at NCBI...but i have the half of all my read (several millons) so i will verify if they mapp on a "transcriptom sequences" just to see if i have a huge part of degraded mRNA in my reads

[westerman]

samtools view fichier bam.bam | awk '$3 =="*"' | samtools view -bS - > no_maps.bam

Hard to explain the above because the 'samtools view fichier' part does not make any sense. Getting rid of the 'fichier' part then means:

1) You are viewing the file 'bam.bam' and converting it from binary BAM format to human readable SAM format.
2) From the SAM choose all lines with the 3rd field equal to "*" (which means no mapping).
3) Convert the SAM file back to BAM and save it.

However a more simple method would be:

samtools view -b -f 4 > no_maps.bam

'-b' means output in BAM format. '-f 4' means to select all reads marked as 'the query sequence itself is unmapped'.

and after i will convert this BAM file into Fastq file with :

java -jar -Xmx6g ~/bin/SamToFastq.jar INPUT=no_maps.bam F=no_map.fastq VALIDATION_STRINGENCY=LENIENT

Well I usually put the '-jar' command next the the jar file itself. Basically to me the '-jar' and '-Xmx6g' parameters are switched around in my mind but I do not think it makes a difference. Basically you are running a Java program that will do the conversion for you. That is just the way Java programs are run.

It looks like you are running a Picard tools program in which case the more modern way of doing this is:

 PicardCommandLine SamToFastq INPUT=no_maps.bam F=no_map.fastq VALIDATION_STRINGENCY=LENIENT

[/QUOTE]

[westerman]

@GenoMax's solution (just capture the unmapped reads first) is good. Redoing the mapping might take some time and you won't learn about samtools nor picard that way but it is what I would have done as the first step. Of course I would have suggested using BBMap instead of bowtie. Am surprised GenoMax didn't suggest that. :-)

[elitatoun]

@ westerman
thank you also for your reply !!!

i think i'm folling in love of this forum !!!

it's clear know for me

i will come back with the results i hope tomorow....
and i will try also your solution...

Why BBmap instead bowtie ?

[GenoMax]

Of course I would have suggested using BBMap instead of bowtie. Am surprised GenoMax didn't suggest that. :-)

Just going with the flow of the original question

"fichier" is "file" so I assume that is just from some instructions @elitatoun was using.

@Elitatoun: This is small non-coding RNA data so perhaps depending on how the prep was done you may see an impact on the % alignment you are going to get.

[elitatoun]

ok thank you Genomax

[elitatoun]

Dear all,
thank you for your advices !
I have finished all my file treatments...and by your advice all it works...

i try to know if the unmapped reads come from degradation of long RNA and particularly mRNA....
so i have mapped the unmapped reads to reference sequence of transcriptome...but i'm not sur about this file (the transcriptome file)...because i don't know if inside the term of "transcriptome" we are talking about only mRNA or all RNA...
do you know a file with only sequences from mRNA ??

[Brian Bushnell]

If you are unsure whether your transcriptome contains the transcripts of interest, I suggest you align to the genome instead (using a splice-aware aligner). Of course, you could just try mapping the unmapped reads to the genome; if they don't map, your transcriptome should be OK.

Unlike a genome, the contents of a transcriptome are not well defined; the accuracy and completeness varies a lot. I think they generally do not contain uRNAs, for example.

[GenoMax]

i try to know if the unmapped reads come from degradation of long RNA and particularly mRNA....
so i have mapped the unmapped reads to reference sequence of transcriptome...but i'm not sur about this file (the transcriptome file)...because i don't know if inside the term of "transcriptome" we are talking about only mRNA or all RNA...
do you know a file with only sequences from mRNA ??

What genome are you working on? Reason I am asking is to see if there is a relatively good known transcriptome for it.

Even though you may have had some degradation aligners should still be able to find a match for those reads if they really belong in your genome. Since they did not align the first time around I am suspicious of the strategy of doing a second round alignment with those unmapped reads to the transcriptome. Perhaps you need to try different options for the original aligner to see if you can get a better result.

What was the result of the blast search for those unmapped reads? Are they mapping to the genome you are working with?

[elitatoun]

@GenoMax,

thank you fr your reply.

I am working with the baus taurus genome.
I will try this afternoon the blast.
But you're right it could be a problem of "options" with the first mapper

We have used bowtie (inside miRdeep2) the first time ==> 50% of the reads unmapped
I have used also bwa (just to see the proportion of mapped reads on the ref genome outside of a miR environnement) ==> 10% of unmapped reads.

so i will work with the options of bowtie !

I have obtained "normal" results at the end (comparison to pubmed) in terms of miR number or putatives....but only 12% of all my reads correspond to miR...and may be i can obtained better results if i modify my extraction procedure...that why i want to know if i have a lot of small reads in the same nucleotide range of my miR, but belonging to degraded mRNA !

I will keep you in touch !
Eli

[GenoMax]

How many cycles was this data (bp per read)?

It would be reasonable to see if blast gives reasonable full length "hits" across the entire read.

You may also want to look at BBMap (if you were planning to do additional alignments).