+/*
+ * getReadsInSnps:
+ * This program takes a set of snps in a custom tab format, and a set of short mapped reads, and evaluates
+ * the sequencing overlap over those snps. Additionally, a miaxture model is fit and used to classify the
+ * snps as homozygous or heterozygous.
+ *
+ * In the final report, the output is:
+ * <snp id> <chromosome> <position> <reference base> <a count> <c count> <g count> <t count> <total count> <snp call>
+ * where snp call is one of:
+ * -1: no call was made (not enough examples to make a call)
+ * 0: the snp is homozygous
+ * 1: the snp is heterozygous
+ *
+ */
#include <sys/types.h>
-#include <dirent.h>
-#include <errno.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <gsl/gsl_statistics.h>
#define WINDOW 25
+#define PI 3.14159265358979323846
+
+//#define DEBUG
using namespace std;
void split (const string& text, const string& separators, vector<string>& words);
char *strrevcomp(const string& input);
+double norm_prob(double x, double mu, double s) { return (1.0)/(s*sqrt(2*PI)) * exp(-0.5*(x-mu)*(x-mu)/(s*s)); }
+
class Loci {
public:
string chr;
};
+
class Read : public Loci {
public:
string seq;
Read& operator=(const Read& r) { this->chr = r.chr; this->pos = r.pos; this->seq = r.seq; return *this;}
};
+typedef vector<Read> Reads;
+
class SNP : public Loci {
public:
string name;
char reference_base;
+ char consensus[4]; // represent the consensus sequence in order. Most often, only the first 1 or 2 will matter.
unsigned int A;
unsigned int C;
unsigned int G;
SNP(string name, string chr, unsigned int pos, char reference_base) : Loci(chr,pos) {
this->name = name;
- this->A = 0; this->C = 0; this->G = 0; this->T = 0; this->total = 0;
+ this->A = 0;
+ this->C = 0;
+ this->G = 0;
+ this->T = 0;
+ this->N = 0;
+
this->reference_base = reference_base;
}
return *this;
}
+ void eval_consensus() {
+ // if A is the max
+ if(A >= C & A >= G & A >= T) { consensus[0] = 'A';
+ if(C >= G & C >= T) { consensus[1] = 'C';
+ if(G >= T) { consensus[2] = 'G'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'G'; }
+ } else if(G >= C & G >= T) { consensus[1] = 'G';
+ if(C >= T) { consensus[2] = 'C'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'C'; }
+ } else { consensus[1] = 'T';
+ if(C >= G) { consensus[2] = 'C'; consensus[3] = 'G'; }
+ else { consensus[2] = 'G'; consensus[3] = 'C'; }
+ }
+
+
+ // if C is the max
+ } else if(C >= A & C >= G & C >= T) { consensus[0] = 'C';
+ if(A >= G & A >= T) { consensus[1] = 'A';
+ if(G >= T) { consensus[2] = 'G'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'G'; }
+ } else if(G >= A & G >= T) { consensus[1] = 'G';
+ if(A >= T) { consensus[2] = 'A'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'A'; }
+ } else { consensus[1] = 'T';
+ if(A >= G) { consensus[2] = 'A'; consensus[3] = 'G'; }
+ else { consensus[2] = 'G'; consensus[3] = 'A'; }
+ }
+ } else if(G >= A & G >= C & G >= T) { consensus[0] = 'G';
+ if(A >= C & A >= T) { consensus[1] = 'A';
+ if(C >= T) { consensus[2] = 'C'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'C'; }
+ } else if(C >= A & C >= T) { consensus[1] = 'C';
+ if(A >= T) { consensus[2] = 'A'; consensus[3] = 'T'; }
+ else { consensus[2] = 'T'; consensus[3] = 'A'; }
+ } else { consensus[1] = 'T';
+ if(A >= C) { consensus[2] = 'A'; consensus[3] = 'C'; }
+ else { consensus[2] = 'C'; consensus[3] = 'A'; }
+ }
+ } else { consensus[0] = 'T';
+ if(A >= C & A >= G) { consensus[1] = 'A';
+ if(C >= G) { consensus[2] = 'C'; consensus[3] = 'G'; }
+ else { consensus[2] = 'G'; consensus[3] = 'C'; }
+ } else if(C >= A & C >= G) { consensus[1] = 'C';
+ if(A >= G) { consensus[2] = 'A'; consensus[3] = 'G'; }
+ else { consensus[2] = 'G'; consensus[3] = 'A'; }
+ } else { consensus[1] = 'G';
+ if(A >= C) { consensus[2] = 'A'; consensus[3] = 'C'; }
+ else { consensus[2] = 'C'; consensus[3] = 'A'; }
+ }
+ }
+ }
+
void add_read(char nuc) {
switch(nuc) {
case 'a':
}
total++;
}
+
+ void clean(unsigned int threshold) {
+ if(A <= threshold) { A = 0; }
+ if(C <= threshold) { C = 0; }
+ if(G <= threshold) { G = 0; }
+ if(T <= threshold) { T = 0; }
+ total = A + C + G + T;
+ eval_consensus();
+ }
+
+ double RE(unsigned int th = 2) {
+ if(total == 0) { return 0.0; }
+
+ double pA = (double)( ((A<th)?A:0)+1e-10)/(double)total;
+ double pC = (double)( ((C<th)?C:0)+1e-10)/(double)total;
+ double pG = (double)( ((G<th)?G:0)+1e-10)/(double)total;
+ double pT = (double)( ((T<th)?T:0)+1e-10)/(double)total;
+
+ //assume equal distribution of A,C,G,T
+ double l2 = log(2);
+ return pA*log(pA/0.25)/l2 + pC*log(pC/0.25)/l2 + pG*log(pG/0.25)/l2 + pT*log(pT/0.25)/l2;
+ }
+};
+
+typedef vector<SNP> SNPs;
+
+//Class to calulate mixture model. Very not general right now, but should be easy enough to make more general
+//if the need arises
+class GaussianMixture {
+
+public:
+ double p;
+ double u1;
+ double s1;
+ double u2;
+ double s2;
+ double Q;
+
+ unsigned int N;
+
+ double delta;
+
+ GaussianMixture(SNPs& snps, double delta = 1e-10) {
+ //initialize model
+ this->p = 0.5;
+ //model 1: heterozygous
+ this->u1 = 1.0;
+ this->s1 = 0.5;
+
+ //model 2: homozygous
+ this->u2 = 2.0;
+ this->s2 = 0.5;
+
+ this->delta = delta;
+ }
+
+ bool classify(double x) {
+ return(norm_prob(x,u1,s1) >= norm_prob(x,u2,s2)) ;
+ }
+
+ // Use EM to fit gaussian mixture model to discern heterozygous from homozygous snps
+ void fit(SNPs& snps, unsigned int count_th) {
+ //initialize relative entropy and probabilities
+ vector<double> RE;
+ vector<double> pr;
+ for(unsigned int i = 0; i < snps.size(); ++i) {
+ if(snps[i].total >= 8) {
+ RE.push_back(snps[i].RE(count_th));
+ pr.push_back(0.5);
+ }
+ }
+
+ this->N = RE.size();
+
+ cerr << this->N << " snps checked\n";
+
+ //calculate initial expectation
+ this->Q = 0.0;
+ for(unsigned int i = 0; i < N; ++i) {
+ Q += pr[i] * (log( this->p ) - log(sqrt(2.0*PI)) - log(this->s1) - (RE[i] - this->u1)*(RE[i] - this->u1)/(2.0*this->s1*this->s1));
+ Q += (1.0-pr[i]) * (log(1-this->p) - log(sqrt(2.0*PI)) - log(this->s2) - (RE[i] - this->u2)*(RE[i] - this->u2)/(2.0*this->s2*this->s2));
+ }
+
+ cerr << "Q: " << this->Q << endl;
+
+ double Q_new = 0;
+ //expectation maximization to iteratively update pi's and parameters until Q settles down.
+ while(1) {
+ cerr << "loop Q: " << Q << endl;
+ Q_new = 0.0;
+
+ double p_sum = 0.0, q_sum = 0.0, u1_sum = 0.0, u2_sum = 0.0;
+ for(unsigned int i = 0; i < N; ++i) {
+ pr[i] = pr[i]*norm_prob(RE[i],this->u1,this->s1) /
+ (pr[i]*norm_prob(RE[i],this->u1,this->s1) + (1.0 - pr[i])*(norm_prob(RE[i],this->u2,this->s2)));
+
+ p_sum += pr[i];
+ q_sum += (1.0 - pr[i]);
+
+ u1_sum += pr[i]*RE[i];
+ u2_sum += (1.0 - pr[i])*RE[i];
+
+ Q_new += pr[i] * (log( this->p ) - log(sqrt(2*PI)) - log(this->s1) - (RE[i] - this->u1)*(RE[i] - this->u1)/(2.0*this->s1*this->s1));
+ Q_new += (1.0-pr[i])* (log(1-this->p) - log(sqrt(2*PI)) - log(this->s2) - (RE[i] - this->u2)*(RE[i] - this->u2)/(2.0*this->s2*this->s2));
+ }
+
+ //update variables of the distributions (interwoven with pi loop to save cpu)
+ this->p = p_sum / this->N;
+ this->u1 = u1_sum / p_sum;
+ this->u2 = u2_sum / q_sum;
+
+ double s1_sum = 0.0, s2_sum = 0.0;
+ for(unsigned int i = 0; i < N; ++i) {
+ s1_sum += pr[i] * (RE[i] - this->u1)*(RE[i] - this->u1);
+ s2_sum += (1.0-pr[i]) * (RE[i] - this->u2)*(RE[i] - this->u2);
+ }
+
+ this->s1 = sqrt(s1_sum/p_sum);
+ this->s2 = sqrt(s2_sum/q_sum);
+
+ if(fabs(this->Q - Q_new) < 1e-5) { break; }
+ this->Q = Q_new;
+ }
+ cerr << "Q: " << Q << endl;
+ }
+
+ void print_model() {
+ cout << "Q: " << Q << " p: " << p << " norm(" << u1 << "," << s1 << ");norm(" << u2 << "," << s2 << ")" << endl;
+ }
};
+
ostream &operator<<( ostream &out, const SNP &h ) {
out << h.name.c_str() << "\t" << h.chr.c_str() << "\t" << h.pos << "\t" << h.reference_base << "\t" << h.A << "\t" << h.C << "\t" << h.G << "\t" << h.T << "\t" << h.total;
+
return out;
}
-typedef vector<Read> Reads;
-typedef vector<SNP> SNPs;
void read_snps(const char* filename, SNPs& snps) {
string delim("\t");
cerr << "Found and sorted " << snps.size() << " snps." << endl;
}
+
+
void read_align_file(char* filename, Reads& features) {
string delim(" \n");
string location_delim(":");
bool strand = ((fields[4].c_str())[0] == 'F')?0:1;
string seq;
- if(strand == 0) {
- seq = fields[0];
- } else {
- seq = strrevcomp(fields[0]);
- }
+ if(strand == 0) { seq = fields[0]; } else { seq = strrevcomp(fields[0]); }
Read read(chr,pos,seq);
features.push_back(read);
}
}
int main(int argc, char** argv) {
- if(argc != 3) { cerr << "Usage: " << argv[0] << " snp_file read_file\n"; exit(1); }
+ if(argc != 4) { cerr << "Usage: " << argv[0] << " snp_file read_file non_reference_output_file\n"; exit(1); }
char snp_filename[1024]; strcpy(snp_filename,argv[1]);
char read_filename[1024]; strcpy(read_filename,argv[2]);
+ char nonref_filename[1024]; strcpy(nonref_filename,argv[3]);
SNPs snps; read_snps(snp_filename, snps);
Reads reads; read_align_file(read_filename, reads);
count_read_at_snps(snps, reads);
+ //fix a guassian mixture model to the snps to classify
+ GaussianMixture g(snps);
+ g.fit(snps, 2);
+
+#ifdef DEBUG
+ g.print_model();
+#endif
+
+ ofstream nonref(nonref_filename);
+ int group;
for(SNPs::iterator i = snps.begin(); i != snps.end(); ++i) {
- cout << (*i) << endl;
+ group = -1;
+ if(i->total >= 10) { i->eval_consensus(); group = g.classify(i->RE()); }
+ cout << (*i) << "\t" << group << "\t";
+ if(group == 0) cout << i->consensus[0];
+ else if(group == 1) cout << i->consensus[0] << "," << i->consensus[1];
+
+ if( ( group == 0 && i->consensus[0] != toupper(i->reference_base) ) || group == 1) {
+ //detected difference from consensus sequence
+ nonref <<i->chr << "\t" << i->pos << "\t";
+ if(group == 0) { nonref << i->consensus[0] << endl; }
+ if(group == 1) {
+ if(i->consensus[0] != toupper(i->reference_base)) {
+ nonref << i->consensus[0] << endl;
+ } else {
+ nonref << i->consensus[1] << endl;
+ }
+ }
+ }
+ cout << endl;
}
+ nonref.close();
}
void split (const string& text, const string& separators, vector<string>& words) {
projects / htsworkflow.git / blobdiff