@shenwei356 when I use this
join df1.fa.tsv df2.fa.tsv | awk '{gsub($3, tolower($3)); print $1"\t"$2}'
it does not work, it seems like there must be both data having the same number of accession , other wise it won't work
I am trying to find a region in a sequence . lets say I have this sequence
>sp|P14306|CPYI_YEAST Carboxypeptidase Y inhibitor OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) OX=559292 GN=TFS1 PE=1 SV=2
MNQAIDFAQASIDSYKKHGILEDVIHDTSFQPSGILAVEYSSSAPVAMGNTLPTEKARSK
PQFQFTFNKQMQKSVPQANAYVPQDDDLFTLVMTDPDAPSKTDHKWSEFCHLVECDLKLL
NEATHETSGATEFFASEFNTKGSNTLIEYMGPAPPKGSGPHRYVFLLYKQPKGVDSSKFS
KIKDRPNWGYGTPATGVGKWAKENNLQLVASNFFYAETK
I want to find the amino acid number of NTLIEYMGPAPPKGSG and make this region lower case
Or remove the amino acid of this 144 to 159 from the sequence
I made a better example
lets call this one df1.txt
>tr|A0A1B1L9R9|A0A1B1L9R9_BACTU ABC transporter permease OS=Bacillus thuringiensis OX=1428 GN=berB PE=4 SV=1
MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVG
MESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIA
ITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLNKELFLKINIVGGLIFLVVSAYSFF
FSCICNDERKALSYSASLTILFFVLDMVGKLSDKLEWMKNLSLFTLFRPKEIAEGAYNIW
PVSIGLIAGALCIFIVAIVVFKKRDLPL
>sp|O15304|SIVA_HUMAN Apoptosis regulatory protein Siva OS=Homo sapiens OX=9606 GN=SIVA1 PE=1 SV=2
MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDE
GCAVVHLPESPKPGPTGAPRAARGQMLIGPDGRLIRSLGQASEADPSGVASIACSSCVRA
VDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
I have a second text file like below lets call it df2
>sp|O15304|SIVA_HUMAN (87-93) Best selected with ratio of 10 , 12 and 14
IGPDGR
>tr|A0A1B1L9R9|A0A1B1L9R9_BACTU (135-168) Not selected with ratio of (12,14)
NKELFLKINIVGGLIFLVVSAYSFF
FSCICNDERKALSYSASLTILFFVLDMVGKLSDKLEWM
I want to merge them both in one file. So few keys are into merging them The output, I am trying to make is like below
>tr|A0A1B1L9R9|A0A1B1L9R9_BACTU ABC transporter permease OS=Bacillus thuringiensis OX=1428 GN=berB PE=4 SV=1 (135-168) Not selected with ratio of (12,14)
MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVG
MESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIA
ITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLnkelflkinivgggliflvvsaysf
ffscicnderkalsysasltilffvldmvgklewmKNLSLFTLFRPKEIAEGAYNIW
PVSIGLIAGALCIFIVAIVVFKKRDLPL
>sp|O15304|SIVA_HUMAN Apoptosis regulatory protein Siva OS=Homo sapiens OX=9606 GN=SIVA1 PE=1 SV=2 (87-93) Best selected with ratio of 10 , 12 and 14
MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDE
GCAVVHLPESPKPGPTGAPRAARGQMLigpdgrLIRSLGQASEADPSGVASIACSSCVRA
VDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
So, first we look at the df1.txt, we find a name like |A0A1B1L9R9| we match the data in the df2.txt and it has the same name. The we find the region that is specific in df2 in parentheses like (87-93) and make it lower case
One solution using SeqKit and shell commands join and awk.
Converting fasta files to tab-delimited format. Only FASTA ID instead of whole header is kept. Sorting is needed if you want join using shell command join later.
$ cat df1.fa | seqkit sort | seqkit fx2tab -i > df1.fa.tsv
$ cat df2.fa | seqkit sort | seqkit fx2tab -i > df2.fa.tsv
$ cat df1.fa.tsv
sp|O15304|SIVA_HUMAN MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDEGCAVVHLPESPKPGPTGAPRAARGQMLIGPDGRLIRSLGQASEADPSGVASIACSSCVRAVDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
tr|A0A1B1L9R9|A0A1B1L9R9_BACTU MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVGMESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIAITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLNKELFLKINIVGGLIFLVVSAYSFFFSCICNDERKALSYSASLTILFFVLDMVGKLSDKLEWMKNLSLFTLFRPKEIAEGAYNIWPVSIGLIAGALCIFIVAIVVFKKRDLPL
Joining two TSV files. Using csvtk if IDs of records in the two files are not equal.
# or:
# csvtk join -H -t df1.fa.tsv df2.fa.tsv
$ join df1.fa.tsv df2.fa.tsv
sp|O15304|SIVA_HUMAN MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDEGCAVVHLPESPKPGPTGAPRAARGQMLIGPDGRLIRSLGQASEADPSGVASIACSSCVRAVDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET IGPDGR
tr|A0A1B1L9R9|A0A1B1L9R9_BACTU MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVGMESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIAITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLNKELFLKINIVGGLIFLVVSAYSFFFSCICNDERKALSYSASLTILFFVLDMVGKLSDKLEWMKNLSLFTLFRPKEIAEGAYNIWPVSIGLIAGALCIFIVAIVVFKKRDLPL NKELFLKINIVGGLIFLVVSAYSFFFSCICNDERKALSYSASLTILFFVLDMVGKLSDKLEWM
Converting subsequence (3rd column) found in sequence (2nd column) to lowercase, and removing 3rd column.
$ join df1.fa.tsv df2.fa.tsv \
| awk '{ gsub($3, tolower($3)); print $1"\t"$2}'
sp|O15304|SIVA_HUMAN MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDEGCAVVHLPESPKPGPTGAPRAARGQMLigpdgrLIRSLGQASEADPSGVASIACSSCVRAVDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
tr|A0A1B1L9R9|A0A1B1L9R9_BACTU MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVGMESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIAITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLnkelflkinivggliflvvsaysfffscicnderkalsysasltilffvldmvgklsdklewmKNLSLFTLFRPKEIAEGAYNIWPVSIGLIAGALCIFIVAIVVFKKRDLPL
Converting TSV back to FASTA format.
$ join df1.fa.tsv df2.fa.tsv \
| awk '{ gsub($3, tolower($3)); print $1"\t"$2}' \
| seqkit tab2fx \
> df1.mask.fa
$ cat df1.mask.fa
>sp|O15304|SIVA_HUMAN
MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDE
GCAVVHLPESPKPGPTGAPRAARGQMLigpdgrLIRSLGQASEADPSGVASIACSSCVRA
VDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
>tr|A0A1B1L9R9|A0A1B1L9R9_BACTU
MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVG
MESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIA
ITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLnkelflkinivggliflvvsaysff
fscicnderkalsysasltilffvldmvgklsdklewmKNLSLFTLFRPKEIAEGAYNIW
PVSIGLIAGALCIFIVAIVVFKKRDLPL
(Optional) If you want keep the full header.
$ seqkit seq -n df1.fa | sed 's/ /\t/' > id2head.tsv
$ cat df1.mask.fa | seqkit replace -k id2head.tsv -p '(.+)' -r '$1 {kv}'
>sp|O15304|SIVA_HUMAN Apoptosis regulatory protein Siva OS=Homo sapiens OX=9606 GN=SIVA1 PE=1 SV=2
MPKRSCPFADVAPLQLKVRVSQRELSRGVCAERYSQEVFEKTKRLLFLGAQAYLDHVWDE
GCAVVHLPESPKPGPTGAPRAARGQMLigpdgrLIRSLGQASEADPSGVASIACSSCVRA
VDGKAVCGQCERALCGQCVRTCWGCGSVACTLCGLVDCSDMYEKVLCTSCAMFET
>tr|A0A1B1L9R9|A0A1B1L9R9_BACTU ABC transporter permease OS=Bacillus thuringiensis OX=1428 GN=berB PE=4 SV=1
MNKQLFLASLKETQKSILSYACGAALYLWLLIWIFPSMVSAKGLNELIAAMPDSVKKIVG
MESPIQNVMDFLAGEYYSLLFIIILTIFCVTVATHLIARHVDKGAMAYLLATPVSRVQIA
ITQATVLILGLLIIVSVTYVAGLVGAEWFLQDNNLnkelflkinivggliflvvsaysff
fscicnderkalsysasltilffvldmvgklsdklewmKNLSLFTLFRPKEIAEGAYNIW
PVSIGLIAGALCIFIVAIVVFKKRDLPL
@shenwei356 when I use this
join df1.fa.tsv df2.fa.tsv | awk '{gsub($3, tolower($3)); print $1"\t"$2}'
it does not work, it seems like there must be both data having the same number of accession , other wise it won't work
Step 2. Joining two TSV files. Using csvtk if IDs of records in the two files are not equal
@shenwei356 I am getting this error [INFO] read key-value file: id2head.tsv [ERRO] no valid data in key-value file: id2head.tsv
checking is id2head.tsv blank.
if it is, paste result of
seqkit seq -n df1.fa | head -n 5 | cat -A
if not, paste result of
head -n 5 id2head.tsv | cat -A
Given sequences in in.fa, e.g.:
>sp|P14306|CPYI_YEAST Carboxypeptidase Y inhibitor OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) OX=559292 GN=TFS1 PE=1 SV=2
MNQAIDFAQASIDSYKKHGILEDVIHDTSFQPSGILAVEYSSSAPVAMGNTLPTEKARSK
PQFQFTFNKQMQKSVPQANAYVPQDDDLFTLVMTDPDAPSKTDHKWSEFCHLVECDLKLL
NEATHETSGATEFFASEFNTKGSNTLIEYMGPAPPKGSGPHRYVFLLYKQPKGVDSSKFS
KIKDRPNWGYGTPATGVGKWAKENNLQLVASNFFYAETK
You could pipe them to the following script lc_pattern.py:
#!/usr/bin/env python
import sys
import re
try:
pattern = sys.argv[1]
except ValueError as ve:
sys.stderr.write("Error: Specify pattern\n")
sys.exit(-1)
def process_record(header, needle):
sequence = needle
for idx in [m.start() for m in re.finditer(pattern, needle, re.IGNORECASE)]:
offset = idx + len(pattern)
sequence = sequence[:idx] + needle[idx:offset].lower() + sequence[offset + 1:]
sys.stdout.write('%s\n%s\n' % (header, sequence))
sequence = ''
for line in sys.stdin:
element = line.rstrip()
if element.startswith('>'):
header = element
if len(sequence) > 0:
process_record(header, sequence)
sequence = ''
else:
sequence += element
# process last record
process_record(header, sequence)
For example:
$ python ./lc_pattern.py NTLIEYMGPAPPKGSG < in.fa > out.fa
$ more out.fa
>sp|P14306|CPYI_YEAST Carboxypeptidase Y inhibitor OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) OX=559292 GN=TFS1 PE=1 SV=2
MNQAIDFAQASIDSYKKHGILEDVIHDTSFQPSGILAVEYSSSAPVAMGNTLPTEKARSKPQFQFTFNKQMQKSVPQANAYVPQD
DDLFTLVMTDPDAPSKTDHKWSEFCHLVECDLKLLNEATHETSGATEFFASEFNTKGSntlieymgpappkgsgHRYVFLLYKQP
KGVDSSKFSKIKDRPNWGYGTPATGVGKWAKENNLQLVASNFFYAETK
@Alex Reynolds
what if I have 100 of them in one file ? is there a possibility to do that . I made an example in my question
This script will process a pattern in a 100 sequences in one FASTA file.
If you have multiple patterns, see the modified script in my answer.
So you could then specify multiple patterns in a pipeline, like so:
$ ./lc_pattern.py ${PATTERN_1} < in.fa | ./lc_pattern.py ${PATTERN_2} | ./lc_pattern.py ${PATTERN_3} | ... | ./lc_pattern.py ${PATTERN_N} > out.fa
I think using re.IGNORECASE should take care of overlapping patterns.
@Alex Reynolds the problem is that I have two txt files and the patterns from one sequence should be matched to the same exact sequences .... This way, I should basically extract all patterns and then do them one by one?
Is there a way that I extract all those seq and do it consecutively without pasting each time?
Log in to answer this question.
@Alex Reynolds I do appreciate your help very much. I am just trying to figure out how to solve it. it has taken so much time from me. I really did not mean to offend you by any mean. You Rock. I appreciate your help for every single word you say. I wish I was able to communicate with you through email. Do you think it is possible?
I followed your post from the start. First, you didn't describe your problem properly, and at each solution Alex Reynolds provided, you added more requirements. This is terrible, because it wastes both your and his time. Be thoughtful when asking questions, lay out the problem in detail from the start, and provide an example of the input data and intended output. A great resource is Tutorial: How To Ask Good Questions On Technical And Scientific Forums.
If you do, upvote his answers / comments, when useful. It is the best way to show appreciation. And his answer and comments were useful - had you asked a good question from the beginning, I bet he would solve your problem right away.
You should use shenwei's answer. I've never been a huge fan of frameworks to handle text files, but his kit is good stuff.
@Alex Reynolds ok thanks
You're very welcome!