Hello,
I wish to filter from a big FASTA file only sequences related to chrs 1-22, X and Y. An example of FASTA sequence is:
>ENSG00000119314|ENST00000210227|PTBP3|9|-1|115024785
GGGTGGCAGGTGCCTGTAATCCCAGCTACTCCAGAGGCTGAGGCAGGGGAATTGCTTGAG
CCTGGGAGGCAGAGGTTGCAGTGAGCCGAGATTGTGCCACTGCACTCCAGCCTGGAGTCT
CACTTTGTCACACAGGGTGGAGTGCAGTGGTGTGATCTCGGCTCACTGCAACCTCTGCTT
ACCGGGTTGAGATTCTCCTGTCTCAACCTCCTGAGTAGCTGGGATTACAGGCGTGCACCA
CCAAGCCAGACTAATTTTCCTATTTTTAGTAGAGATGGGGGGTTTCACCATGTTGGCCAG...
>ENSG00000236011|ENST00000211377|GPANK1|HSCHR6_MHC_COX|-1|31616421
CCCTATTCCTACCTAACCTCCCCTCAGGACTCAGGCTCCAATGTGTTGAGCCCCAACTCC
TTCCCATAAGACTGCCACACGGTGCTTTCCTTTCCCTTCTTCAACACTCACCAATGGGAA
GCATTGGCTGGTTCTCACAGTACACACGAGGACAGTAACCAAAGTCTCCTTGCTGGTACT
TTTCCAACTGAGGTGAATACAATGGAAGGGGTTGGCAGGTAGATGTAAAGAAGAGGCAAC
TCCCTTCGCAGCCCAACCCATACCACTCTGTCCCCCACTCCTCCCACCTCTGTCCAGAGG
CCCCTTCTCTGGACTAGACGGGCTCTCAAACTTCTGTGTTGCCTTTCTTCCAATTAGGCA
GGCTACAAACCATCAGAGCCATTTGTTGTTTGTTCCTTGAGGAAGAGGCAGTCTATCACA
ACTCTCTGATTCAAGGTCTGTCTCCCTCCCTGAAAACAATCCCTTCAGGATGACCCCCAA...
>ENSG00000087494|ENST00000201015|PTHLH|12|-1|28115255
TCCGCTCACGGGCCCCGAGACCCCCGAAGTTCCCATGGAGCCTAAGATCCCCAGGAGCCA
AGCCTGCCCCGTCCCTGCGGATCAGCTTCCTAATGGGCGACCCAAGTCTATCGCAGGCGG
TGGGGATGAGGACGCTGGGTGGGAGGAGGGGAGGGGAGGCTGAAAAAGATCATCCCCCTT
GCCCTAAGGCCTCTCCCAAGACCCTGGACCCCTGCCCTAAGAGACTCAGGCCTCCCTTGC
TGCAGTGGGAGCGCAAACACCAGGGCAGGAGACTCCAGAGAAGGAGCGCATAACTCAACG
TTTGCTCTCCTGAAGCCTTATTTCTGATAAAAATTACAGAAAAGTTAGGCAGGATCCAAA
GACACCGTAATGACCAGCTCAAAGCCAAACAGACAGGACATCCAGTGCGGGTGTCTGGAT...
As you can see in the forth place after "|" there is the chromosome name: in the first one is: 9 the second one is: HSCHR6MHCCOX and the third one is: 12
I want to create a new FASTA file containing only 9, 12 sequences and in more generally 1-22, X and Y sequences.
What is fastest way to do it?
Thanks,
Tom.
•
written
by
sthait •
5 answers
just one quick observation, the correct way to write the segment of the chromosome is chr2:1-2000, not chr2:1:2000
THank you, leszek. It is the most convenient and comparatively fastest way I have found to remove the scaffolds in the reference genome.
You can do the following:
chr=$(seq 1 1 22)" X Y"
for i in $chr ; do
grep ">" input.fa | grep "|$i|" | awk 'BEGIN {FS=">"} {print $2}' > selection_$i
samtools faidx input.fa `cat selection_$i` > output_$i.fa
done
This is a bash script, so you can just run it as is in a terminal.
It did not work for me - nothing was printed to file. However this modification did:
for i in $chr ; do
samtools faidx $fasta $i >> chr_sequences.fa
done
It should be noted that this puts all sequences into a single file.
What is $fasta here? You did not make this clear...
$fasta is a variable that represents the input fasta. It should be defined before the loop as fasta=my_fancy.fasta
if you want Ensemble gene models, you can use BioMart to obtain your sequence filtered: http://uswest.ensembl.org/biomart/martview
You can do some pattern matching using pyfaidx:
faidx --regex "^.+\|.+\|.+\|[[1-9][0-2]?|[XY]]" input.fa > output.fa
•
written
by
Matt Shirley
Log in to answer this question.