Question: Is there a way to automate or pipline bedtofasta or similar script?
0
gravatar for alexjfortuna
2.6 years ago by
alexjfortuna0 wrote:

Hi,

This is my first post so please bear with me. I have about 2000 bed files that I would like to extract FASTA sequences from. I have been able to extract a single set of bed sequences using the following command:

.narrowPeak is a BED file.

bedtools getfasta -fi TAIR10_Chr.fasta -bed /copy_peaks/ABI5_col_v3h.narrowPeak > ABI5_col_v3h.fasta

within my copy_peaks folder there are 2000 other files I would like to perform this operation on. Any ideas?

The results will then be used in fimo, where again I will need to perform this operation on the extracted fasta files.

Thanks,

sequencing assembly genome • 1.1k views
ADD COMMENTlink modified 2.6 years ago by WouterDeCoster40k • written 2.6 years ago by alexjfortuna0
3
gravatar for Alex Reynolds
2.6 years ago by
Alex Reynolds28k
Seattle, WA USA
Alex Reynolds28k wrote:

I wrote a Perl script a while back that uses samtools to automate this conversion, using indexed FASTA files:

If you use bash, you could write a for loop to automate conversion of a batch of files.

ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by Alex Reynolds28k
3
gravatar for WouterDeCoster
2.6 years ago by
Belgium
WouterDeCoster40k wrote:

I think a bash loop could help you, or for more advanced options you could use gnu-parallel.

If all your files have the same .narrowPeak extension you could use the following example:

for f in /copy_peaks/*.narrowPeak
do
r=`echo $f | cut -f 3 -d '/' | sed s/narrowPeak/fasta/`
bedtools getfasta -fi TAIR10_Chr.fasta -bed $f > $r
done
ADD COMMENTlink modified 2.6 years ago • written 2.6 years ago by WouterDeCoster40k

I'll give this a shot, all of files have the same extension. Does this also create a new file? and I was hoping for the new file to keep the name of the narrowpeak file but with .fasta as the suffix. Thanks!

ADD REPLYlink written 2.6 years ago by alexjfortuna0

The first thing the loop does is create a $r variable, which is based on the file name of the input as you can see. That $r variable is then used as output name and indeed I switch the narrowPeak to fasta and write to the current directory. Changing this should be pretty straightforward, right?

You can test the output names first using the following loop:

for f in /copy_peaks/*.narrowPeak
do
r=`echo $f | cut -f 3 -d '/' | sed s/narrowPeak/fasta/`
echo $r
done
ADD REPLYlink written 2.6 years ago by WouterDeCoster40k

I ran this :

for f in /copy_peaks/*.narrowPeak
do
r=`echo $f | cut -f 3 -d '/' | sed s/narrowPeak/fasta/`
echo $r
done

The result was: *.fasta I then attempted to run the loop as you had written it, the results were:

for f in /copy_peaks/*.narrowPeak; do r=`echo $f | cut -f 3 -d '/' | sed s/narrowPeak/fasta/`; bedtools getfasta -fi /mnt/c/Users/Desktop/FIMO_Project/TAIR10_seq_chr/TAIR10_Chr.fasta -bed $f > $r; done

bash: *.fasta: No such file or directory

Please pardon all my questions, I've been stuck on this for a while now.

ADD REPLYlink modified 2.6 years ago by WouterDeCoster40k • written 2.6 years ago by alexjfortuna0

I assume based on this that there is no path /copy_peaks/. Are you sure of that path? Is it perhaps just copy_peaks in the current directory, and not on the root directory?

ADD REPLYlink written 2.6 years ago by WouterDeCoster40k

Yes, copy_peaks is a folder in the current working directory, now I'm getting $r is an ambiguous redirect? Thanks - Alex

ADD REPLYlink written 2.6 years ago by alexjfortuna0

Does this look okay? Note that I changed the first line to reflect the directory in which the bed files can be found.

for f in copy_peaks/*.narrowPeak
do
r=`echo $f | cut -f 2 -d '/' | sed s/narrowPeak/fasta/`
echo $r
done

Your example bedtools getfasta -fi TAIR10_Chr.fasta -bed /copy_peaks/ABI5_col_v3h.narrowPeak > ABI5_col_v3h.fasta was then incorrect because the path to the narrowPeak file doesn't exist as such. Should have been bedtools getfasta -fi TAIR10_Chr.fasta -bed copy_peaks/ABI5_col_v3h.narrowPeak > ABI5_col_v3h.fasta

ADD REPLYlink modified 2.6 years ago • written 2.6 years ago by WouterDeCoster40k

It's worked! Thank you for your help

ADD REPLYlink written 2.6 years ago by alexjfortuna0
Please log in to add an answer.

Help
Access

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