Question: How to extract a peptide sequence (interval) from a FASTA file?
1
gravatar for pedro19ms
12 days ago by
pedro19ms20
pedro19ms20 wrote:

Hello,

I have two files:

  1. sequences.fasta
  2. intervals.txt

I used the "seqinR" package to obtain the protein sequences of interest (~ 3500 proteins sequences) and saved them to "sequences.fasta":

> sp | O94827 | PKHG5_HUMAN
MDDQSPAEKKGLRCQNPACMDKGRAAKVCHHADCQQLHRRGPLNLCEACDSKFHSTMHYDGHVRFDLPPQG
SVLARNVSTRSCPPRTSPAVDLEEEEEESSVDGKGDRKSTGLKLSKKKARRRHTDDPSKECFTLKFDLNVDIETEIVPAMKKKSLGEVLLP
VFERKGIALGKVDIYLDQSNTPLSLTFEAYRFGGHYLRVKAPAKPGDEGKVEQGMKDSKSLSLPILRPAGTGPPAL
ERVDAQSRRESLDILAPGRRRKNMSEFLGEASIPGQEPPTPSSCSLPSGSSGSTNTGDSWKNRAASRFSGFFSS
GPSTSAFGREVDKMEQLEGKLHTYSLFGLPRLPRGLRF

> sp | P10515 | ODP2_HUMAN
MWRVCARRAQNVAPWAGLEARWTALQEVPGTPRVTSRSGPAPARRNSVTTGYGGVRALCGWTPSSGATPRNRLLLQLL
GSPGRRYYSLPPHQKVPLPSLSPTMQAGTIARWEKKEGDKINEGDLIAEVETDKATVGFESLEECYMAKILVAEGTRDVPIGA
IICITVGKPEDIEAFKNYTLDSSAAPTPQAAPAPTPAATASPPTPSAQ

My "intervals.txt" file has the protein ID, the starting and ending position of my peptide of interest:

ID               start         end
O94827           1             69
P10515           2             120

I would like to know if there is any way to obtain in a third .txt file the protein ID and the sequence that corresponds to that interval? I would appreciate it if someone helped me!

Thanks in advance!

R • 124 views
ADD COMMENTlink modified 12 days ago by bioinformatics2020340 • written 12 days ago by pedro19ms20

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.
code_formatting

ADD REPLYlink written 12 days ago by RamRS28k
2
gravatar for bioinformatics2020
12 days ago by
bioinformatics2020340 wrote:
library(Biostrings)
fa <- readBStringSet("file.fa")
intervals <- read_tsv("intervals.txt")
for(i in seq_along(fa)){ 
  fa[[i]] <- subseq(fa[[i]], start = intervals$start[[i]], end = intervals$end[[i]])
}
writeXStringSet(fa, "file.fa")
ADD COMMENTlink written 12 days ago by bioinformatics2020340

Thank very you! I made some minor adjustments and it worked perfectly!

ADD REPLYlink written 11 days ago by pedro19ms20
2
gravatar for cpad0112
12 days ago by
cpad011213k
India
cpad011213k wrote:

with seqkit:

input:

$ cat test.fa
> sp | O94827 | PKHG5_HUMAN
MDDQSPAEKKGLRCQNPACMDKGRAAKVCHHADCQQLHRRGPLNLCEACDSKFHSTMHYDGHVRFDLPPQG
SVLARNVSTRSCPPRTSPAVDLEEEEEESSVDGKGDRKSTGLKLSKKKARRRHTDDPSKECFTLKFDLNVDIETEIVPAMKKKSLGEVLLP
VFERKGIALGKVDIYLDQSNTPLSLTFEAYRFGGHYLRVKAPAKPGDEGKVEQGMKDSKSLSLPILRPAGTGPPAL
ERVDAQSRRESLDILAPGRRRKNMSEFLGEASIPGQEPPTPSSCSLPSGSSGSTNTGDSWKNRAASRFSGFFSS
GPSTSAFGREVDKMEQLEGKLHTYSLFGLPRLPRGLRF

> sp | P10515 | ODP2_HUMAN
MWRVCARRAQNVAPWAGLEARWTALQEVPGTPRVTSRSGPAPARRNSVTTGYGGVRALCGWTPSSGATPRNRLLLQLL
GSPGRRYYSLPPHQKVPLPSLSPTMQAGTIARWEKKEGDKINEGDLIAEVETDKATVGFESLEECYMAKILVAEGTRDVPIGA
IICITVGKPEDIEAFKNYTLDSSAAPTPQAAPAPTPAATASPPTPSAQ


$ cat test.txt 
O94827  1   69
P10515  2   120

output:

$ awk  -v OFS="\t" 'NR>1 {print $1,$2-1,$3}' test.txt > test.bed
$ seqkit seq -w 0 test.fa | seqkit subseq --quiet --bed test.bed - --id-regexp "\|\s([A-z0-9]+)\s\|" -w 0
>O94827_1-69:. 
MDDQSPAEKKGLRCQNPACMDKGRAAKVCHHADCQQLHRRGPLNLCEACDSKFHSTMHYDGHVRFDLPP
>P10515_2-120:. 
WRVCARRAQNVAPWAGLEARWTALQEVPGTPRVTSRSGPAPARRNSVTTGYGGVRALCGWTPSSGATPRNRLLLQLLGSPGRRYYSLPPHQKVPLPSLSPTMQAGTIARWEKKEGDKIN
ADD COMMENTlink modified 11 days ago • written 12 days ago by cpad011213k
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: 1562 users visited in the last hour