fastq.cpp

#include <iostream>
#include <stdlib.h>
#include <fstream>
#include <string.h>
#include <sstream> // for int2strinng
#include <vector>
#include <algorithm> // for count
#include <numeric> // for accumulate
#include <map> // for sorted map.
#include <tr1/unordered_map>
using namespace std;
#define DEFAULT -999

const char *usage = "============================================\n"
    "fastq  : utilities for summarizing and filtering fastq files\n"
    "author : brentp (bpederse@gmail.com)\n"
    "license: MIT\n"
    "\n"
    "usage:\n\t$ fastq [action] [options]\n"
    "where [action] is one of 'summarize' or 'filter':\n"
    "\t$ fastq summarize [summarize_options]\n"
    "      [summarize_options] are:\n"
    "\t--codon : print a codon usage table with columns A,C,T,G,N and one\n"
    "\t\trow per read basepair\n"
    "\t--quality : print a quality table with columns:\n"
    "\t\tmean, std, median, 1, 2.5, 5, 10 , 25, 75, 90, 95, 97.5, 99%\n"
    "\t\tand one row per read-basepair\n"
    "\n"
    "     [filter_options] are\n:"
    "\t--quality [X] mask bases with quality < X to 'N'\n"
    "\t--n [X] exclude records with more than X 'N's (after above masking)\n"
    "\t\tdefault 1000 is no masking. 4 is a reasonable value\n"
    "\t--unique after above filtering, print only unique sequences.\n"
    "\t(for those groups with more than one record, choose the one\n"
    "\t with the best quality\n. NOTE: order is not maintained\n."
    "\n"
    "except when called with --unique\n"
    "both 'summarize' and 'filter' require an additional option:\n"
    "\t--adjust [N] where N is likely either 33 or 64. see:\n"
    "\thttp://en.wikipedia.org/wiki/FASTQ_format\n"
;
void error(const char *msg){
    fprintf(stderr, "**%s\n", msg);
}

typedef struct fastq_record {
    string seq;
    string header;
    string squality;
    string oheader;
    vector<int> iquality;
    float average_quality;
    int ns;
    bool operator< (const fastq_record &b) const{
        return average_quality < b.average_quality;
    };
} fastq_record;

inline void write_fastq_record(fastq_record fq){
    printf("%s\n%s\n%s\n%s\n", 
                    fq.header.c_str(),
                    fq.seq.c_str(),
                    fq.oheader.c_str(),
                    fq.squality.c_str());
}

typedef struct opts {
    string action;
    int adjust;
    bool filter_unique;
    int filter_quality;
    int filter_n;
    bool summarize_codon;
    bool summarize_quality;
    string filename;
} opts;

typedef vector<fastq_record> record_list;

// http://www.cplusplus.com/forum/articles/13355/
// http://www.codeguru.com/cpp/cpp/cpp_mfc/stl/article.php/c15303__2
//typedef std::tr1::unordered_multimap<string, fastq_record> record_multimap;
typedef tr1::unordered_map<string, record_list> record_map;
typedef tr1::unordered_map<char, int> char_int_map;
typedef tr1::unordered_map<int, char_int_map> bp_char_freq_map;

typedef tr1::unordered_map<int, int> int_int_map;
typedef tr1::unordered_map<int, int_int_map> bp_qual_freq_map;

void add_quality(fastq_record *qr, int adjust, int filter_limit, int *masked){
    int l = qr->squality.length();
    qr->iquality.resize(l);
    qr->squality.resize(l);

    for(int i=0; i < l; i++){
        qr->iquality[i] = qr->squality[i] - adjust;
        if(qr->iquality[i] < filter_limit){
            qr->seq[i] = 'N';
            *masked += 1;
        }
    }
    qr->average_quality = accumulate(qr->iquality.begin(), 
                                     qr->iquality.end(), 0) / (float)l;
    
}

string int2string(int number) {
   stringstream ss; //create a stringstream
   ss << number; //add number to the stream
   return ss.str(); //return a string with the contents of the stream
}

int chart_adjust(int nreads){
    return nreads > 10000000 ? 10000 : nreads > 1000000 ? 1000 : nreads > 100000 ? 100 : nreads > 10000 ? 10 : 1;
}


int fsummarize_quality(opts o){
    bp_qual_freq_map quality_map;
    int_int_map imap;
    struct fastq_record qr;
    ifstream fh (o.filename.c_str(), ios::in);
    if (!fh.is_open()){ 
        fprintf(stderr, "does '%s' exist? not able to open\n", 
                o.filename.c_str());
        return 1;
    }
    float _stat_pts[] = {0.01, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.99};
    vector<float> stat_pts (_stat_pts, _stat_pts + sizeof(_stat_pts) / sizeof(float));
    vector<float>::iterator stat_it;

    int i, L, nreads=0, masked=0;
    while (1){
        getline(fh, qr.header);
        if(fh.eof()){ break; }
        getline(fh, qr.seq);
        getline(fh, qr.oheader);
        getline(fh, qr.squality);
        add_quality(&qr, o.adjust, 0, &masked);
        L = qr.iquality.size();
        for(i=0; i < L; ++i){
            // so here we have a lookup of read position -> quality_score -> count.
            quality_map[i][qr.iquality[i]]++; 
        }
        ++nreads;

    }
    long qsum_quality;
    int qi;

    vector<int> counts;
    string url = "http://chart.apis.google.com/chart?chbh=a&chdlp=l&cht=lc&chs=800x300&chco=FF0000,00FF00,0000FF,FFFF00,CCCCCC";
    //url += "&chxt=x,x&chxr=0,1," + int2string(l) + "," + int2string(l - 1) + "&chxl=1:|read position|&chxp=1," + int2string(2 * l / 3);
    url += "&chd=t0:";
    for(i=0; i < L; ++i){
        qsum_quality = 0;
        imap = quality_map[i];
        counts.clear();
        for(int_int_map::const_iterator it=imap.begin(); it != imap.end(); ++it){
            // first is the quality. 2nd is the count.
            qsum_quality += (it->first * it->second);
            // faster way to do this?
            for(qi=0; qi < it->second; qi++){ counts.push_back(it->first); }
        }

        // sort then find the pctiles.
        sort(counts.begin(), counts.end());

        string stats = "";
        for(stat_it=stat_pts.begin(); stat_it != stat_pts.end(); stat_it++){
            stats += int2string(counts[nreads * (*stat_it)]) + "\t"; 
        }

        printf("%i\t%s%.5f\n", i + 1, stats.c_str(), (float)qsum_quality / nreads);
    }
    return 0;
}

int fsummarize_codon(opts o){
    bp_char_freq_map cm;
    struct fastq_record qr;
    ifstream fh (o.filename.c_str(), ios::in);
    if (!fh.is_open()){ 
        fprintf(stderr, "does '%s' exist? not able to open\n", 
                o.filename.c_str());
        return 1;
    }
    // this is o.summarize_codon
    int i, l, si;
    int nreads = 0;
    while (1){
        getline(fh, qr.header);
        if(fh.eof()){ break; }
        getline(fh, qr.seq);
        getline(fh, qr.oheader);
        getline(fh, qr.squality);
        l = qr.squality.length();
        for(i=0; i < l; ++i){
            // default is 0, so no need to intialize.
            // hash of read position -> char -> count.
            cm[i][qr.seq[i]]++;
        }
        nreads += 1;
    }
    string actgn = "ACGTN";
    string url = "http://chart.apis.google.com/chart?chbh=a&chdl=A|C|G|T|N&chdlp=l&cht=bvs&chs=800x300&chco=FF0000,00FF00,0000FF,FFFF00,CCCCCC";
    url += "&chxt=x,x&chxr=0,1," + int2string(l) + "," + int2string(l - 1) + "&chxl=1:|read position|&chxp=1," + int2string(2 * l / 3);
    int div = chart_adjust(nreads);
    url += "&chds=0," + int2string(nreads/ div + 1) + "&chd=t:";
    printf("#reads: %i\n", nreads);
    for(si=0; si < actgn.size(); si++){
        for (i=0; i < l; ++i){
            int val = cm[i][actgn[si]] / div;
            url += int2string(val) + (i == l - 1 ? "" : ",");
        }
        if(si < actgn.size() -1) {url += "|";}
    }

    printf("#%s\n", url.c_str());
    printf("read_pos\tA\tC\tG\tT\tN\n");
    for (i=0; i < l; ++i){
        printf("%i", i + 1);
        for(si =0; si < actgn.size(); si++){
            int val = cm[i][actgn[si]];
            printf("\t%i", val);
        }
        printf("\n");
    }
    return 1;
}

int write_filtered(record_map &rm, opts &o){
    fprintf(stderr, "in write_filtered\n");
    record_list rl;
    int i;

    for(record_map::const_iterator it = rm.begin(); it != rm.end(); ++it){
        // it->first is the sequence key.
        rl = it->second;
        // only write the one with best avg_q
        if(o.filter_unique){
            // sort by average_quality.
            sort(rl.rbegin(), rl.rend());            
            write_fastq_record(rl[0]);
        }            
        else {
            for(i=0; i < rl.size(); i++){
                write_fastq_record(rl[i]);
            }
        }
    } 
}


int filter(opts o){
    record_map rm;
    struct fastq_record qr;
    int ns_skipped = 0;
    int i = 0;
    int quality_skipped = 0;
    ifstream fh (o.filename.c_str(), ios::in);
    if (!fh.is_open()){ 
        cout << "does '" << o.filename << "' exist? not able to open" << endl;
        return 1;
    }
    record_map::iterator pos;
    while (1){
        getline(fh, qr.header);
        if(fh.eof()){ break; }
        i += 1;
        getline(fh, qr.seq);
        getline(fh, qr.oheader);
        getline(fh, qr.squality);

        add_quality(&qr, o.adjust, o.filter_quality, &quality_skipped);

        qr.ns = count(qr.seq.begin(), qr.seq.end(), 'N');
        if (qr.ns > o.filter_n){
            ns_skipped += 1;
            continue;
        }
        
        pos = rm.find(qr.seq);
        if (pos == rm.end()){
            // create a new vector.
            vector<fastq_record> rl;
            rl.push_back(qr);
            rm.insert(record_map::value_type(qr.seq, rl));
        }
        else {
            // already seen, add it to the vector.
            rm[qr.seq].push_back(qr);
        }
    }
    fprintf(stderr, "masked %i basepairs because of quality's\n", 
            quality_skipped);
    float ratio = (float)ns_skipped/i * 100.0;
    fprintf(stderr, "skipped %i reads because of N's(%.3f %%)\n", 
            ns_skipped, ratio);

    ratio = (float)rm.size()/i * 100.0;
    fprintf(stderr, "read %i unique of %i total records(%.3f %%)\n",
            (int)rm.size(), i, ratio);
    return write_filtered(rm, o);
}



int parse_args(int argc, char *argv[], opts *o){
    int i;
    o->filter_n = 1000;
    o->filter_unique = false;
    o->filter_quality = DEFAULT;
    o->summarize_codon = false;
    o->summarize_quality = false;
    o->adjust = DEFAULT;
    o->action = (string)argv[1];
    o->filename = (string)argv[argc - 1];
    
    if(o->action == "summarize"){
        for(i=2; i < argc - 1;){
            if(string(argv[i]) == "--codon"){
                o->summarize_codon = true;
            }
            else if (string(argv[i]) == "--quality"){
                o->summarize_quality = true;
            }
            else if (string(argv[i]) == "--adjust"){
                o->adjust = atoi(argv[i + 1]);
                i += 1;
            }
            else {
                fprintf(stderr, "bad argument: %s\n", argv[i]);
            }
            i++; 
        }
    }
    //bool filter_unique;
    //int filter_quality;
    else if (o->action == "filter"){
        for(i=2; i < argc - 1;){
            if(string(argv[i]) == "--unique"){
                o->filter_unique = true;
            }
            else if (string(argv[i]) == "--adjust"){
                o->adjust = atoi(argv[i + 1]);
                i += 1;
            }
            else if (string(argv[i]) == "--quality"){
                o->filter_quality = atoi(argv[i + 1]);
                i += 1;
            }
            else if (string(argv[i]) == "--n"){
                o->filter_n = atoi(argv[i + 1]);
                i += 1;
            }
            else {
                fprintf(stderr, "bad argument: %s\n", argv[i]);
            }
            i++;
        }

    }
    else {
        return 1;
    }
    if(o->adjust == -999 && o->filter_quality != DEFAULT && o->action == "filter"){
        error("must send in value for --adjust");
        return 1;
    }
    return 0;
}

int main(int argc, char *argv[]){
    if (argc < 3){
        error(usage);
        return 1;
    }
    opts o;
    if(parse_args(argc, argv, &o) != 0){
        error(usage);
        return 1;
    }

    record_map rm;

    if (o.action == "filter"){
        filter(o);
    }
    else if (o.action == "summarize"){
        if (o.summarize_codon){ fsummarize_codon(o); }
        if (o.summarize_quality) { fsummarize_quality(o); }
    }

    return 0;
}