Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • lh3
    replied
    @dpryan your reheader is better in the long run. Haven't thought about that...

    Leave a comment:


  • dpryan
    replied
    If you have a BAM file, the chromosome/contig names are only stored in the header (there's an ordered list and each read just has a numeric value saying where in that list it's chromosome name can be found). Yeah, if you have a bunch of different mappings, then following what Heng Li suggested might prove faster.

    Leave a comment:


  • rkizen
    replied
    Thanks dpryan!

    I have not used reheader before. Will it replace all the sequence names? or will it just replace the header info?

    There are a lot of contigs and there are many mappings which may not have the same contigs represented. I will have to write a script to run "samtools view -H file.bam > header.sam" for each mapping and then replace the sequence names for each of the header sam files and apply those back onto the origin bam file using samtools reheader.

    For now I just need some intervals quickly so I will use the method lh3 to get these intervals and then work on converting the headers later.

    Leave a comment:


  • dpryan
    replied
    The reheader script can simply be "samtools reheader"

    Do a "samtools view -H file.bam > header.sam", edit the header.sam file by replacing the chromosome names (they need to be in the same order as the original, so don't reorder anything!), then use samtools reheader. Unless your alignments are to a large number of contigs, this should be quite quick.

    Leave a comment:


  • rkizen
    replied
    Thanks lh3! That works. I should change the sequence names but I needed this for something a bit urgent and there are a lot of sequence names which would require a bit of time ensuring that a name changing script was not corrupting the file in any way. Thank you again and I am kicking myself for not using the quotation marks to get past the pipes in the sequence names.

    Leave a comment:


  • lh3
    replied
    Try "gi|358485511|ref|NC_006088.3|:100-1000", though I don't know if that works.

    The best solution is to replace the sequence names with something easier.

    Leave a comment:


  • rkizen
    started a topic Samtools view bam parse region issues

    Samtools view bam parse region issues

    I have mapped reads from illumina sequencing onto a reference genome using bwa aln and I wanted to get some regions of the bam files using samtools view:

    samtools view -o test_region.sam merged.sort.bam 1:100-1000

    however I kept getting the error:

    [bam_parse_region] fail to determine the sequence name.

    I have realized that I used an NCBI version of the reference genome which resulted in the reference sequence name in the bam file to be that of the NCBI sequence names. For example, here is bit of the bam file:

    EBRI093151_0050:4:117:5648:11807#0 99 gi|358485511|ref|NC_006088.3| 66454 60 101M = 66940 587 CTCCAAGATCATCCAGTCCACCCATCCACCCACCACCAATGTAACCCCATTAAAACACGTCCCTCAGTACCACATCTAAATGTTTCTTGAGCACTGGCAGG fffffdffffefffffefeffffeffffffffffffffcdefedffdfeffdfaeffffdeffffdedfffeffe`f`ccdfadfffdaYd^cd`bbd`dc XT:A:U NM:i:0 SM:i:37 AM:i:37 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:101


    Notice how instead of "chr1" it has "gi|358485511|ref|NC_006088.3|" as the sequence name. I believe this is the cause of the parsing issue. However, now I cannot seem to figure out how to pull a region from the bam file. I tried using:

    samtools view -o test_region.sam merged.sort.bam NC_006088.3:1-1000

    But I get the same error as before:
    [bam_parse_region] fail to determine the sequence name.

    If I put this instead:

    samtools view -o test_region.sam merged.sort.bam NC_006088.3

    It seems to do as one would expect and pull all read mappings for that sequence region, however, I would like to pull specific coordinates from each chromosome.

    Does anybody have some ideas of how I could do this? I was thinking about replacing all the reference sequence names with their generic equivalent (ie "chr1,chr2 etc") however I am not certain that would solve my issue here.

    Thank you!

Latest Articles

Collapse

  • seqadmin
    Best Practices for Single-Cell Sequencing Analysis
    by seqadmin



    While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
    06-06-2024, 07:15 AM
  • seqadmin
    Latest Developments in Precision Medicine
    by seqadmin



    Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

    Somatic Genomics
    “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
    05-24-2024, 01:16 PM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 06:54 AM
0 responses
10 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-14-2024, 07:24 AM
0 responses
15 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-13-2024, 08:58 AM
0 responses
14 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-12-2024, 02:20 PM
0 responses
17 views
0 likes
Last Post seqadmin  
Working...
X