How to extract a peptide sequence (interval) from a FASTA file?

Question

[Peter]

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!

[Ram]

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.

Answer 1

[bioinformatics2020]

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")

[Peter]

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

Answer 2

[cpad0112]

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