SEQanswers

SEQanswers (http://seqanswers.com/forums/index.php)
-   General (http://seqanswers.com/forums/forumdisplay.php?f=16)
-   -   how convert several bam files to fasta? (http://seqanswers.com/forums/showthread.php?t=90651)

Manuelly 08-27-2019 11:31 AM

how convert several bam files to fasta?
 
Hello!
I have 1979 mtDNA bam files and I need to convert these files to individual fastas. Since there are so many files, I would like to do this conversion faster or even once. I have no idea how to do that.
Can someone help me?

GenoMax 08-27-2019 03:52 PM

Use reformat.sh from BBMap suite. Take a look at BAM processing options in in-line help to decide if you want to keep primary reads etc.

Code:

for i in *.bam; do reformat.sh in=${i} out=${i}.fa add_options_you_need; done
Above would get you individual reads.

If your "convert bam to fasta" need is for consensus fasta sequence then check this tutorial on Biostars.

Manuelly 08-28-2019 07:11 AM

so using reformat.sh command can I convert multiple bam to fasta at once? for this i need to open the folder with all bam and run the command?
Sorry, if the question is too obvious, but I'm still very lay.
I have difficulty understanding the command options in the terminal, is there any online manual that I help with BBmap?

GenoMax 08-28-2019 08:30 AM

Yes you can convert multiple files at once. If your data came from paired-end reads then be aware that the "out=" files will contain data in interleaved format. You will need to separate R1/R2 reads. Again reformat.sh can be used for that as well.

Here is a guide for reformat.sh.

Manuelly 09-11-2019 07:44 AM

I tried to execute the command and the following message appeared. Do you know how I can solve this?
/media/lucasmalzoni/Seagate Backup Plus Drive1/Joao/MT/02-Bam$ reformat.sh in=${i}.bam out=${i}.fa
java -ea -Xms300m -cp /home/lucasmalzoni/anaconda3/opt/bbmap-38.67-0/current/ jgi.ReformatReads in=.bam out=.fa
Executing jgi.ReformatReads [in=.bam, out=.fa]

Exception in thread "main" java.lang.RuntimeException: Can't read file '.bam'
at shared.Tools.testInputFiles(Tools.java:1157)
at jgi.ReformatReads.<init>(ReformatReads.java:337)
at jgi.ReformatReads.main(ReformatReads.java:50)
(base) lucasmalzoni@galton:/media/lucasmalzoni/Seagate Backup Plus Drive1/Joao/MT/02-Bam$

SNPsaurus 09-11-2019 10:58 AM

Can you copy what you did as a command and paste it here, and then list the directory as well and paste a few filenames as well? It looks like it didn't find any bam files. Are there files ending in .bam in that directory?

Manuelly 09-12-2019 05:48 AM

1.I opened the directory with the files in the terminal.
2. Then I used the following command:
reformat.sh in=${i}.bam out=${i}.fa
3. I received the following message:
java -ea -Xms300m -cp /home/lucasmalzoni/anaconda3/opt/bbmap-38.67-0/current/ jgi.ReformatReads in=.bam out=.fa
Executing jgi.ReformatReads [in=.bam, out=.fa]

Exception in thread "main" java.lang.RuntimeException: Can't read file '.bam'
at shared.Tools.testInputFiles(Tools.java:1157)
at jgi.ReformatReads.<init>(ReformatReads.java:337)
at jgi.ReformatReads.main(ReformatReads.java:50)
4.file extensions contained in the directory
/home/lucasmalzoni/Documentos/DOC/1-file_aligned_SortSam_MarkDuplicates.bai
/home/lucasmalzoni/Documentos/DOC/1-file_aligned_SortSam_MarkDuplicates.bam
/home/lucasmalzoni/Documentos/DOC/1-file_aligned_SortSam_MarkDuplicates.metrics.txt

I don't understand why it's not working.

SNPsaurus 09-12-2019 08:11 AM

You should do the entire command as Genomax listed (with a change, I think, I'll mention below):

for i in *.bam; do reformat.sh in=${i}.bam out=${i}.fa; done

The first part (for i in *.bam) finds all the files ending in .bam in the directory, then passes them to the "do" command. So you should change the command to the below which removes the .bam from the in= name, since $i should already have the .bam at the end.

for i in *.bam; do reformat.sh in=$i out=${i}.fa; done

The output file will be named something like
filename.bam.fa

GenoMax 09-12-2019 08:23 AM

Good catch SNPsaurus. I have updated my original command. There was an additional ".bam" in there.

@Manuelly: Use
Code:

for i in *.bam; do reformat.sh in=$i out=${i}.fa; done
If you are trying to run this for just one file then you should do:

Code:

reformat.sh in=your.bam out=file.fa

RobertWood90 09-14-2019 11:14 AM

Convert several bam files to fasta? I had no idea about this before as I've only been doing some web development for netticasino-opas.net and kasinobonukset24.net lately. I got full ideas from this article. Looks like I can do it myself. Thank you.

Manuelly 09-17-2019 10:39 AM

thank you very much for all the help. I realize that I have difficulty understanding the syntax of some commands and I don't know how to overcome this difficulty. Like "for i in *.bam". Is there any material I can read and improve my understanding?

Manuelly 09-17-2019 11:50 AM

1.after extracting the consensus in fasta format, I opened fasta and came across this:
>E00382:156:HMFJCCCXX:1:1102:7659:38315/2
TTCCCCTTAAATAAGACATCACGATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGT
>E00382:156:HMFJCCCXX:1:1105:22019:53803/1
ATGGATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCCCA
2.I think my fasta is separated by reads. and to solve this I used the following command:
for i in *.fa; do reformat.sh in1=$i in2=$i out=${i}.fa ; done
3. I believe it is not this command or something went wrong in the command. Can anyone help me?

I want uninterrupted fastas. it is possible?

GenoMax 09-17-2019 12:00 PM

Output fasta reads were interleaved (put in one file) since you had paired-end reads to begin with. To separate them into two files you should use the command below.

Code:

for i in *.fa; do name=$(basename $(i} .fa); reformat.sh in=${name}.fa out1=${name}_R1.fa out=${name}_R2.fa ; done
You could also do the original conversion by providing "out1=" and "out2=" separate file names during BAM conversion to avoid this step.

Manuelly 09-17-2019 12:46 PM

I think I did not express myself well. I need a single fasta per sample and not interleaved. like that:
>1 NC_012920.1:1-16569
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTT
CGTCTGGGGGGTGTGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTC
GCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATT
ACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATA
ACAATTGAATGTCTGCACAGCCGCTTTCCACACAGACATCATAACAAAAAATTTCCACCA
AACCCCCCCTCCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAA
AACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCA
CTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAA

GenoMax 09-17-2019 03:03 PM

Then we have been on wrong track all along. You should have said you need a consensus fasta file from the BAM alignment.

What you need is: https://www.biostars.org/p/367960/

I had posted this in post #2 in this thread above.


All times are GMT -8. The time now is 07:18 PM.

Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.