match two files and filter out

Question

I have two files

file A

>TCONS_00000075 gene=XLOC_000030
CCGCCGGCTGCTGCGCGCACCGACTTGTCACCACCCCAGCACGTCCTCCACGTATACAAGCGCTACGGTC
CACCGCGGCAGCGTCGACGTCCTTGTCCGCAAACATGGTGGTGGCAGCTTCCTCATCGAGCAGCAGCAAC
TCATCCTCGAGGGGAAGGGCCCAGAGCTTCTAATCCTACACGGCAACAACACTTTATACTTGTGTATAAT
>TCONS_00013830 gene=XLOC_006942
AAACACGGTTAGCTTGATATCACTGATGATCGATGGGATAGAGTCAGAGAACATCTTGTTCCTTAATTAT
CTCAATTCGTGAGATGTTGGACGATATCTCGATAGGGAGAGAAGGCGTTGTTCTGGATCATCACCGTGCT
CAGGGGTCAATTTTACACTGAGCAGGGGCAAAGACGTAAATTTTTACTTCCTTACTTGAGTAAGAGCAAG
TTTAATACTACAACCAACTACTACAAACTCCAATTCATTTATAACCAATCTAATAACTTATTCATACAAT
AGTTACCTATAAGCATATACTACACACACAACGTATTGGAATCCTCCGTGCTGCTGCTGGCTACAGATCT

file B

XLOC_000030
XLOC_000059
XLOC_000210

FileA is a FASTA sequence file. Each line starting with > is a sequence name and the lines beneath it are the sequence. I want to extract the sequences of those IDs mentioned in FileB. In this case:

file C

>TCONS_00000075 gene=XLOC_000030
CCGCCGGCTGCTGCGCGCACCGACTTGTCACCACCCCAGCACGTCCTCCACGTATACAAGCGCTACGGTC
CACCGCGGCAGCGTCGACGTCCTTGTCCGCAAACATGGTGGTGGCAGCTTCCTCATCGAGCAGCAGCAAC
TCATCCTCGAGGGGAAGGGCCCAGAGCTTCTAATCCTACACGGCAACAACACTTTATACTTGTGTATAAT

I tried this command:

perl -pe 's/\n//; s/>(.*)/\n>$1\t/' A |grep -f <(awk '{print $1}' B) |sed 's/\t/\n/' | fold -w 60 > C

but it's not working.

Answer 1

#!/bin/bash
while read line
do
        grep -A 1 $line filea >> filec
done < fileb

Answer 2

Tip: Use database tools for database work.

If you end up spending all of your time just working out the mechanics of looking things up as opposed to doing the lookups themselves, and your commands to look things up are long combinations of perl, sed, awk, and grep that start to resemble modem line noise, then it is time to consider using actual database tools for your database instead of cobbling bespoke queries together with text-processing tools.

select * from RNA where gene in ('XLOC_000030','XLOC_000059','XLOC_000210');

is a lot simpler to work with.

That said, have some text processing modem line noise.

Don't duplicate work.

You are transforming your database on the fly from a multiple-line-per-record form into a one-line-per-record form. Do that once, not every query. Have a Makefile that says:

A.flat: A
        sed -e '/^>/s/$$/ /;:a;$$!N;s/\n //;ta;s/^>//;P;D' $^ > $@.tmp
        mv $@.tmp $@

Then just run make every time that A could have changed. (Don't forget that make relies on that whitespace at the start of those lines being a TAB character.)

That sed program is:

# Append a space to the first line in the record.
/^>/s/$$/ /
:a
# Join if not EOF
$!N
# Eliminate newline and space if the joined line began with a space.
s/\n //
# Loop around if so.
ta
# Eliminate the initial start of record character, as that is now the line break.
s/^>//
# Print and remove the whole record from the buffer.  Then start again.
P
D

Querying

You've over-complicated your query command, which is actually just

grep -f B A.flat > C.flat

Converting the query output back to multiple-line records is more line noise:

sed -e 's/^/>/' C.flat|fold -s -w70|sed -e 's/^[^>]/ /' > C

Omitting the intermediate files thus gives:

grep -f B A.flat|sed -e 's/^/>/'|fold -s -w70|sed -e 's/^[^>]/ /' > C

And again …

To repeat: Actual database tools would be much better. If you have a large number of records, which seems likely given what you appear (from this and other questions) to be doing, searches, insertions, and deletions are going to be very inefficient using the text-processing tools method. An actual database could be indexed on the gene field, in contrast.