Extract nucleotides from fasta file by a single position

Question

[misbahabas]

Asslamu Alikum,

I have a tab separated file like

chr38 456

chr38 555

and I want fetch the nucleotide at this position from fasta file. Is there any perl script or linux command to do this work?.
Please tell me if have any idea about it

thankx in advance

[misbahabas]

thanks for reply

Actually I have only single position of SNP not starting and ending position, file not in bed format

[WouterDeCoster]

Please use ADD COMMENT or ADD REPLY to answer to previous reactions, as such this thread remains logically structured and easy to follow. I have now moved your post but as you can see it's not optimal. Adding an answer should only be used for providing a solution to the question asked.

I said pretty close, I guess you can easily make a bed format from your file.

[misbahabas]

Bed file require starting and ending position , but i have only one position , sorry i am new in bioinformatics , i have no idea about it

kindly tell me how make bed file from this single position

[WouterDeCoster]

You are going to enjoy Cheat Sheet For One-Based Vs Zero-Based Coordinate Systems

[cpad0112]

seqtk subseq in.fa reg.bed > out.fa from https://github.com/lh3/seqtk

Answer 1

[WouterDeCoster]

bedtools getfasta comes pretty close to what you ask for

Answer 2

[Matt Shirley]

You can format the table like UCSC region strings (chr:start-end) which are 1-based closed coordinates. This way you don't have to worry about converting to [0,1) half-open coordinates like BED requires. Then you can pipe the region strings to either samtools faidx or faidx from the pyfaidx package.

pip install pyfaidx
awk '{printf "%s:%s-%s\n",$1,$2,$2}' < regions.txt | xargs faidx ref.fasta

Answer 3

[st.ph.n]

A Python solution. Linearize your fasta file, if it has more than one line per sequence. Save this code as get_pos.py. Change the names of pos.txt to your tab-delimited file, and input.fasta to the name of your fasta file. Run as python get_pos.py > nucl_at_pos.txt. nucl_at_pos.txt can be changed to whatever output file name you wish.

#!/usr/bin/env python

with open('pos.txt', 'r') as f:
    pos = {}
    for line in f:
        pos[line.strip().split('\t')[0]] = int(line.strip().split('\t')[1])

with open('input.fasta', 'r') as fh:
    seqs = {}
    for line in fh:
        if line.startswith(">"):
            seqs[line.strip().split('>')[1]] = next(fh).strip()

for i in pos:
    print i, '\t', seqs[i][pos[i]]

[Matt Shirley]

Just a suggestion for simplifying this process: you don't need to linearize the fasta file, and you don't need to read the whole thing into memory - especially if you're just extracting a few single nucleotides.

[st.ph.n]

Thanks for the tips, Matt. I should learn more about pyfaidx. Having a single-line fasta to work with is just a preference of mine, otherwise I would read the file in with Biopython. I like also to write parts out so the OP can learn something and understand what's going on, and without knowing the size of either file, I agree can be overly taxing for both time and computing power when simpler options are available.

This part could also be changed to this, to iterate through the sequences instead of loading to memory, again depending on the size of the files:

with open('input.fasta', 'r') as fh:
    for line in fh:
        if line.startswith(">"):
            chrom = line.strip().split('>')[1]
            if chrom in pos:
                print chrom, '\t', pos[chrom], '\t', next(f).strip()[pos[chrom]]

[Matt Shirley]

It's good to do the file parsing yourself, especially as a learning exercise. I think for production or publication it's nice to use parsers that are tested and hopefully eliminate a broad class of simple bugs.

Answer 4

[drkennetz]

for this line:

 chrom = line.strip().split('>')[1]

It could be replaced with:

chrom = line.strip().strip('>').split()[0]

Which would get rid of the '>' and just return the value (for example when running on an actual genome). printing this line for genome hg19 yields: 1 10 11 12 13 14 15 16 17 18 19 2 20 21 22 3 4 5 6 7 8 9 MT X Y.

They are all printed to new line. The next issue with the parsing above is that if you have a large list of key:value pairs in your pos dictionary you will run into a problem with duplicates. It will not increment keys, something else to think about although I haven't worked on a fix just yet.

Answer 5

[Rashedul Islam]

This should work.

awk '{print $1 "\t" $2-1 "\t" $2}' input.file | bedtools getfasta -fi genome.fa -bed stdin -fo output.fa

bedtools's getfasta is here.