How to extract a peptide sequence (interval) from a FASTA file?
Hello,
I have two files:
sequences.fastaintervals.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!
• 2,059 views
•
link
updated
by
bioinformatics2020
•
written
by
Peter •
•
written
by
Peter •
2 answers
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")
• 0 views
•
link
written
by
bioinformatics2020
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
Log in to answer this question.
Please use the formatting bar (especially the
codeoption) to present your post better. You can use backticks for inline code (`text` becomestext), 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.