Question: (Closed) knowing the exact position where each read is mapped
0
gravatar for A
4.4 years ago by
A3.7k
A3.7k wrote:

hi biostars,

by tophat2 " tophat2 -p 8 -G genes.gtf genome file.fastq" command (gtf i think is the annotation flie and genome is the whole genome fasta), i produced a file named accepted_hits.bam which using "

bam2bed < accepted_hits.bam  | bedmap --echo --count genes.bed - > answer4.bed


" command,  by beddops tool or

bedtools multicov -bams yourBAMfile.bam -bed yourGenes.bed command 
by bed tool,  i have a bed file now, like below

I    334    646    "YAL069W    .    +    protein_coding    CDS    0    exon_number "1"; gene_id "YAL069W"; gene_name "YAL069W"; p_id "P3633"; protein_id "YAL069W"; transcript_id "YAL069W"; transcript_name "YAL069W"; tss_id "TSS1128";    2

for example 2 reads have been mapped on gene YAL069W in 334-646 distance...

but how i can know where each of these reads individually have been mapped???? i mean how i can have a file in which i could see where each read individually has been mapped...for example these two reads in protein_coding (CDS) part, where have been mapped individually not just knowing that in 334 - 646, two reads have been mapped on gene YAL069W (not a distance in where a number of reads have been mapped)....

please if you have any idea, let me know

thank you so much

rna-seq tophat2 ribo-seq • 1.2k views
ADD COMMENTlink modified 4.4 years ago by Annika Forsingdal210 • written 4.4 years ago by A3.7k
1

"""

how i can have a file in which i could see where each read individually has been mapped

"""

You never encountered SAM format in your life? Samtools do not ring the bell? RTFM

ADD REPLYlink written 4.4 years ago by Darked894.2k

hey don't be angry please...if you read my previous posts you could observe that i tried bedops, bedtools, samtools...i have a file with number of reads but this number is the sum of mapped reads in a given distance but i need to know where the each read is mapped individualy.i used samtools

for example "samtools view yourBAMfile.bam chr1:567876-568100 | wc -l" but there is no output just running something in cmd...after running this is the example of output of this command in cmd

XVI    946855    946858    "YPR204C-A    .    -    protein_coding    stop_codon    0    exon_number "1"; gene_id "YPR204C-A"; gene_name "YPR204C-A"; p_id "P5953"; transcript_id "YPR204C-A"; transcript_name "YPR204C-A"; tss_id "TSS5620";    120
XVI    946855    947338    "YPR204C-A    .    -    protein_coding    exon    .exon_number "1"; gene_id "YPR204C-A"; gene_name "YPR204C-A"; p_id "P5953"; seqedit "false"; transcript_id "YPR204C-A"; transcript_name "YPR204C-A"; tss_id "TSS5620";    1927
XVI    946858    947338    "YPR204C-A    .    -    protein_coding    CDS    0exon_number "1"; gene_id "YPR204C-A"; gene_name "YPR204C-A"; p_id "P5953"; protein_id "YPR204C-A"; transcript_id "YPR204C-A"; transcript_name "YPR204C-A"; tss_id "TSS5620";    1922
XVI    947335    947338    "YPR204C-A    .    -    protein_coding    start_codon    0    exon_number "1"; gene_id "YPR204C-A"; gene_name "YPR204C-A"; p_id "P5953"; transcript_id "YPR204C-A"; transcript_name "YPR204C-A"; tss_id "TSS5620";    196
XVI    947698    947701    "YPR204W    .    +    protein_coding    stop_codon    0    exon_number "1"; gene_id "YPR204W"; gene_name "YPR204W"; p_id "P5115"; transcript_id "YPR204W"; transcript_name "YPR204W"; tss_id "TSS840";    7
[izadi@lbox161 bin]$

but again this is just the totall number of the mapped reads in given distance and not clear the exact position of where each read has been mapped

 

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by A3.7k
1

man wc

man less

man head

If you expect that piping samtools view output into wc will give you where the individual reads are mapping then you need to go back and learn what these commands do. Including being able to check that you can execute samtools, that the BAM file is there, it is not corrupted/overwritten by some badly executed command etc. Nobody can help you with this except by dictating line by line and resolving problems ppl have in their first few hrs of working with a command line. Read i.e.:

http://openwetware.org/wiki/Wikiomics:WinterSchool_day1#Introduction_to_Linux_and_the_command_line

And check that you have samtools on the PATH & chr1 in your BAM file header SQ as a sanity check. 

Good luck

ADD REPLYlink modified 4.4 years ago • written 4.4 years ago by Darked894.2k

since i joined biostars never i have faced such an angry biostar!!...im working with linux since about a month ago..anyway thanks

ADD REPLYlink written 4.4 years ago by A3.7k
1

Hello Fereshteh!

We believe that this post does not fit the main topic of this site.

Not a slightest attempt at solving the problem before posting

For this reason we have closed your question. This allows us to keep the site focused on the topics that the community can help with.

If you disagree please tell us why in a reply below, we'll be happy to talk about it.

Cheers!

ADD REPLYlink written 4.4 years ago by Darked894.2k
1
gravatar for Annika Forsingdal
4.4 years ago by
Copenhagen, Denmark
Annika Forsingdal210 wrote:

You can view your bam files in IGV to see to positions of where each read have mapped.

It only requires indexing of your bam file, e.g. using samtools:

$samtools index accepted_hits.bam

ADD COMMENTlink written 4.4 years ago by Annika Forsingdal210
Please log in to add an answer.
The thread is closed. No new answers may be added.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 836 users visited in the last hour