upstream and downstream regions extraction

Question

[harpreetmanku04]

Hello, I have a query related to upstream and downstream regions extraction. See I have done blastn for 30 nt small query sequences against 500 nt long db sequences now my question is after running blastn, how can I extract upstream and downstream regions for my 30 nt small query sequences?

Answer 1

[Sam]

If you have the coordinate of the alignment e.g. loc, then you can extract the sequence from the reference 500nt sequence starting from loc-30 with length 30nt+2*30nt providing that you have no insertion or deletion in the alignment. If you have the indel, then just add the corresponding length.

To actually do the extraction, you can use the substr function from most programme, e.g. R, perl, awk.

[harpreetmanku04]

but i want to extract after the alignment with standalone blast and also how can i come to know that in my local database this is the reference sequence for that small sequence?

[Sam]

After the alignment, you should be able to get the alignment file, which should contain the sequence id for the reference sequence. You can then get the reference sequence fasta file which you can use to extract the required sequence.

[harpreetmanku04]

Yes I got the subject id along with query id. Thank you. And is there any R script or perl script (understandable) for extracting the sequence (particular region)?

[Sam]

I normally do it using awk, but I guess you can try the Biostrings package in R

[harpreetmanku04]

hello, I'm stuck with the blast command, see when I do blastn for 300 bp long seq then it retains result but when I do blastn for 20-25 nucl seq then it retains no result. Do you have any idea why it is so? I want to blast small sequences with my local database.

[Sam]

Maybe you can post a new question to get an answer to that? I am not familiar with blast so I cannot give you definitive answer. However, the most straightforward guess will be that there simply be no sequence similarity with your small sequence and your local database.

[harpreetmanku04]

Can you post awk command?

[harpreetmanku04]

Hello Sam Sir, can you operate my file to extract upstream and downstream regions using your command(s). I shall be thankful to you. actually I tried but I failed and for me it is taking too much time if you could do it or you can provide me awk command I can try that one also. please help me out here.

[Sam]

Assuming your input is a fasta file containing only your target sequence, then you can use this script:

awk -v seq="" '{if(!($1~"^>")){seq=seq$0}else{print seq; seq="";}}END{print substr(seq, loc-30, 30+2*30)}' <input.fa>

Where loc is the coordinate of the alignment. If you need to do the for different sequence and different coordinates, you will need a more complicated script.

The first part of the script is essentially concatenating the whole sequence into one string and the second part (substr) is where you extract the region of interest

[harpreetmanku04]

Actually I have excel file which contains sequence id, start and stop respectively (generated by blastn) position separately and fasta file which contain sequences separately and I want to extract region from that fasta file.

This command gives error.

[harpreetmanku04]

I tried this command but it also generates an error.

awk 'NR==FNR{if($0~/^>/){i=substr($0,2);getline};a[i]=a[i] $0;next}{print ">" $1 ORS substr(a[$1], $2, $3-$2+1)}' file1 FS=\\t file2

[harpreetmanku04]

Can't we connect on team viewer, if you can solve it over there? Thanku :)

[Sam]

If you put your excel content into a text file (tab delimited), you can do the following:

samtools faidx reference.fasta
awk '{print $1">>output.fasta";print "samtools faidx reference.fasta "$2":"$3"-"$4" >> output.fa"}' <input> |bash

Your input file need to be in the format of <id>\t<chr>\t<start>\t<end>

If you start and end position doesn't account for the flanking region, then you can use this

samtools faidx reference.fasta
awk '{print $1">>output.fa"; print "samtools faidx reference.fasta "$2":"$3-30"-"$4+30" >> output.fa"}' <input> | bash

The result will be print to output.fa do delete the file or rename it in the awk code before each run, otherwise the result will only keep appending onto the file

And for your information Extract User Defined Region From An Fasta File

[harpreetmanku04]

do i need to install sam tool for this?

Answer 2

[harpreetmanku04]

no my problem is like this:

I have file1 is this

GL3482.1 GAACTTGAGATCCGGGGA GCAGTGGATCTCCACCAG CGGCCAGAACTGGTGCAC CTCCAGGCCAGCCTCGTC CTGCGTGTC GL3550.1 GGCATTTTTGTGTAATTTTTGGCTGGATGAGGT GACATTTTCATTACTACCATTTTGGAGTACA GL3472.1 TTTTCCTGTTCACTGCTGCTTTTCTATAGACAGCA GCAGCAAGCAGTAAGAGAAAGTA ```

file2 is this:

seq id     start    end
GL3482.1   323100   323743
GL3550.1   41911    40888
GL3472.1   274408   272617

and I want result like this:

>GL3482.1 
TTGAGA
>GL3550.1 
GGCATTTTTGTG
>GL3472.1
TTCCTGTT

[Sam]

samtools faidx reference.fasta
awk '{print ">"$1">>output.fa"; print "samtools faidx reference.fasta "$2":"$3-30"-"$4+30" | tail -n +2 >> output.fa"}' <input> | bash

Your example output actually make no sense to me because:

1. The examples doesn't corresponds to your coordinates (I understand that you want to make the post short).
2. You didn't state if you want the upstream/downstream
3. Do you want to include the actual sequence, or just the up/downstream.

[harpreetmanku04]

Yes I have to include sequence between upstream and downstream region too.

[Sam]

Then what is the problem of the above script?