Here is a bit of code that you can use to extract FASTA files from PLINK genotype calls. I am doing a few assumptions:
- There is a reference genome to compare against
- If no data is available from the genotyping we will infer the reference call. This might be acceptable for some Whole Genome Sequencing data sources (and even so...), but if you are using something less dense then you might want to mark it as missing data. The version without the reference genome would be much easier to code
- Adding to the reference genome, your genes of interest will be on a BED file
- I am doing this inside IPython Notebook. I only make use of a magic, which can be trivially changed by os.system
- Your PLINK is binary encoded
Lets start with the import code, plus some variables that you need to change:
import gzip from Bio.Seq import Seq from Bio import SeqIO from Bio.SeqRecord import SeqRecord indiv_file = 'indivs.txt' genome = '/data/atroparvus/tra/ag/ref/pest3.fa.gz' my_bed = 'my_genes.bed' plink_pref = '/data/gambiae/phase1/AR3'
You will need to change the genome file (a gzipped reference fasta), the bed file with your positions on interest) and the plink prefix.
Next, lets do a simple function to parse BED files, this will be coded as a generator:
def get_features(bed_file): for l in open(bed_file): toks = l.rstrip().split('\t') chrom = toks start = int(toks) end = int(toks) gene = toks strand = toks yield chrom, start, end, gene, strand
Now lets get the sequence for a certain area of interest:
def get_seq(fasta_name, chrom, start, end, determine_chrom=lambda x: x.description.split(':')): recs = SeqIO.parse(gzip.open(fasta_name, 'rt', encoding='utf8'), 'fasta') for rec in recs: if determine_chrom(rec) == chrom: return rec.seq[start:end + 1]
This will take your reference genome file and return a Biopython sequence with your area of interest. Notice the parameter to determine the chromosome from the FASTA description. Different FASTA files can have completely different descriptor lines, so we abstract that away. In my concrete case, the description looks like this:
So, I take the description, split it by : and get the 3rd element (2L in this case). "But I could have used the ID", I hear you say. Yes I could, but I wanted to illustrate code for strange description lines...
OK, now lets get the SNPs from the area of interest in the genome from our PLINK file. The file is in binary format (you can easily change this for text or even VCF as long as you use PLINK 1.9+) and the output is is text format (so that we can parse it). We are using IPython magic here (os.system, if you prefer):
def get_snps(plink_prefix_in, plink_prefix_out, indiv_file, chrom, start, end): !plink --bfile $plink_prefix_in --chr $chrom --from-bp $start --to-bp $end\ --keep $indiv_file --allow-extra-chr --recode --out $plink_prefix_out
OK, now the work-horse of all of this, a function to read a PLINK file and return sequences based on the reference genome:
def get_plink_seqs(ref_seq, plink_pref, start): f = open(plink_pref + '.map') poses =  for l in f: pos = l.rstrip().split('\t')[-1] poses.append(int(pos)) f = open(plink_pref + '.ped') for l in f: toks = l.rstrip().replace('\t', ' ').split(' ') fam = toks ind = toks alleles = toks[6:] seq1 = list(ref_seq) seq2 = list(ref_seq) for i, pos in enumerate(poses): ref_pos = pos - start ref_allele = ref_seq[ref_pos] a1 = alleles[2 * i] a2 = alleles[2 * i + 1] if ref_allele == ref_allele.lower(): #maintain case of the ref a1 = a1.lower() a2 = a2.lower() seq1[ref_pos] = a1 seq2[ref_pos] = a2 #if a1 != ref_allele or a2 != ref_allele: # print(fam, ind, i, ref_pos, pos, ref_allele, a1, a2) yield fam, ind, ''.join(seq1), ''.join(seq2)
So, the input a reference sequence (only the tiny slice of the genome corresponding to your area of interest), the PLINK prefix of your PLINK files and the start position. The code starts by loading the PLINK positions from the MAP file (notice that I ignore the chromosome: you should only have your chromosome, plus PLINK can do strange things to chromosome names - I will not detail that here). Most of the code involves offsetting from the absolute coordinates in the PLINK file to the relative coordinates in our reference sequence. Also note that we maintain the capitalization as in the reference genome (which normally has a meaning, like soft-masking).
Finally, lets extract the sequences proper to FASTA files: these will be file named by the feature/gene:
for feature in get_features(my_bed): chrom, start, end, gene, strand = feature ref_seq = get_seq(genome, chrom, start, end) print(chrom, start + 1, end + 1, gene) get_snps(plink_pref + '/' + chrom, gene, indiv_file, chrom, start + 1, end + 1)# 1 based my_records =  my_records.append(SeqRecord(seq=ref_seq, id=gene, description='Reference')) for fam, ind, seq1, seq2 in get_plink_seqs(ref_seq, gene, start + 1): my_records.append(SeqRecord(seq=Seq(seq1), id=gene, description='%s_%s_1' % (fam, ind))) my_records.append(SeqRecord(seq=Seq(seq2), id=gene, description='%s_%s_2' % (fam, ind))) SeqIO.write(my_records, '%s.fa' % gene, 'fasta')
We iterate through all our features on the BED file: we get the corresponding reference sequence, get the SNPs from PLINK and generate the FASTA sequences. We use SeqIO for writing.
- The PLINK files are not phased, so you cannot assume correct phasing of the FASTA files
- I am assuming diploid data here
- The code could be made more efficient by getting the correct reference sequence only once (it can be quite heavy). This can be easily achieved by retrieving the BED features by chromosome order (many times that have that order, anyway) and only loading the reference sequence when the chromosome/scafold changes<
I hope this code is useful to you. As usual comments and criticism are very welcome.