How to extract a peptide sequence (interval) from a FASTA file?
2
1
Entering edit mode
10 months ago
Peter ▴ 20

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 • 357 views
ADD COMMENT
0
Entering edit mode

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 REPLY
2
Entering edit mode
10 months ago
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 COMMENT
0
Entering edit mode

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

ADD REPLY
2
Entering edit mode
10 months ago

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 COMMENT

Login before adding your answer.

Traffic: 2553 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6