Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • some help ! BAM to Fastq

    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

  • #2
    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.

    Comment


    • #3
      @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

      Comment


      • #4
        Originally posted by elitatoun View Post

        Code:
        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:
        Code:
        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 :

        Code:
        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:

        Code:
         PicardCommandLine SamToFastq INPUT=no_maps.bam F=no_map.fastq VALIDATION_STRINGENCY=LENIENT
        [/QUOTE]

        Comment


        • #5
          @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. :-)

          Comment


          • #6
            @ 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 ?

            Comment


            • #7
              Originally posted by westerman View Post
              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.

              Comment


              • #8
                ok thank you Genomax

                Comment


                • #9
                  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 ??

                  Comment


                  • #10
                    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.

                    Comment


                    • #11
                      Originally posted by elitatoun View Post
                      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?

                      Comment


                      • #12
                        @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

                        Comment


                        • #13
                          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).

                          Comment

                          Latest Articles

                          Collapse

                          • seqadmin
                            Exploring the Dynamics of the Tumor Microenvironment
                            by seqadmin




                            The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
                            07-08-2024, 03:19 PM
                          • seqadmin
                            Exploring Human Diversity Through Large-Scale Omics
                            by seqadmin


                            In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
                            06-25-2024, 06:43 AM

                          ad_right_rmr

                          Collapse

                          News

                          Collapse

                          Topics Statistics Last Post
                          Started by seqadmin, 07-10-2024, 07:30 AM
                          0 responses
                          25 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, 07-03-2024, 09:45 AM
                          0 responses
                          201 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, 07-03-2024, 08:54 AM
                          0 responses
                          211 views
                          0 likes
                          Last Post seqadmin  
                          Started by seqadmin, 07-02-2024, 03:00 PM
                          0 responses
                          193 views
                          0 likes
                          Last Post seqadmin  
                          Working...
                          X