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
"""
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
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 cmdXVI 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
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
since i joined biostars never i have faced such an angry biostar!!...im working with linux since about a month ago..anyway thanks
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!