← Back to branch summary

~ubuntu-branches/ubuntu/utopic/r-cran-genabel/utopic

« back to all changes in this revision

Viewing changes to src/export_plink.cpp

  • Committer: Package Import Robot
  • Author(s): Andreas Tille, Charles Plessy, Andreas Tille
  • Date: 2014-08-07 17:30:04 UTC
  • mfrom: (1.1.6)
  • Revision ID: package-import@ubuntu.com-20140807173004-24c20tmrw61yl5dz
Tags: 1.8-0-1
[ Charles Plessy ]
* debian/control: removed myself from Uploaders.

[ Andreas Tille ]
* New upstream version
* Moved debian/upstream to debian/upstream/metadata
* cme fix dpkg-control
* More detailed copyright

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/export_plink.cpp
5
5
extern "C" {
6
6
#endif
7
7
 
8
 
std::string* getGenotype(std::string coding,std::string sep)
9
 
{
10
 
        std::string* Genotype = new (std::nothrow) std::string [4];
11
 
        std::string Letter0 = coding.substr(0,1);
12
 
        std::string Letter1 = coding.substr(1,1);
13
 
        Genotype[0] = "0"+sep+"0";
14
 
        Genotype[1] = Letter0+sep+Letter0;
15
 
        Genotype[2] = Letter0+sep+Letter1;
16
 
        Genotype[3] = Letter1+sep+Letter1;
17
 
        return Genotype;
18
 
}
19
 
 
20
 
SEXP export_plink(SEXP Ids, SEXP Snpdata, SEXP Nsnps, SEXP NidsTotal, SEXP Coding, SEXP From, SEXP To,
21
 
                SEXP Male, SEXP Traits, SEXP Pedfilename, SEXP Plink, SEXP Append)
22
 
{
23
 
 
24
 
        int from = INTEGER(From)[0];
25
 
        int to = INTEGER(To)[0];
26
 
        std::vector<unsigned short int> sex;
27
 
        sex.clear();
28
 
        unsigned short int sx;
29
 
        for(unsigned int i=(from - 1);i<to;i++) {
30
 
                sx = INTEGER(Male)[i];
31
 
                if (sx==0) sx=2;
32
 
                //Rprintf("%d %d\n",i,sx);
33
 
                sex.push_back(sx);
34
 
        }
35
 
        std::vector<std::string> ids;
36
 
        for(unsigned int i=0;i<((unsigned int) length(Ids));i++)
37
 
                ids.push_back(CHAR(STRING_ELT(Ids,i)));
38
 
 
39
 
        std::vector<std::string> coding;
40
 
        for(unsigned int i=0;i<((unsigned int) length(Coding));i++)
41
 
                coding.push_back(CHAR(STRING_ELT(Coding,i)));
42
 
 
43
 
        //Rprintf("0\n");
44
 
        unsigned int nsnps = INTEGER(Nsnps)[0];
45
 
        int nids = to - from + 1;
46
 
        int nidsTotal = INTEGER(NidsTotal)[0];
47
 
        int ntraits = INTEGER(Traits)[0];
48
 
        bool append = LOGICAL(Append)[0];
49
 
        bool plink = LOGICAL(Plink)[0];
50
 
        std::string filename = CHAR(STRING_ELT(Pedfilename,0));
51
 
        std::ofstream fileWoA;
52
 
        int ieq1 = 1;
53
 
        char * snpdata = (char *) RAW(Snpdata);
54
 
 
55
 
        //      int gtint[nidsTotal];
56
 
        int *gtint = new (std::nothrow) int[nidsTotal];
57
 
 
58
 
        //Rprintf("nsnps=%d\n",nsnps);
59
 
        //Rprintf("nids=%d\n",nids);
60
 
        //Rprintf("to=%d\n", to);
61
 
        //Rprintf("from=%d\n", from);
62
 
 
63
 
        //char gtMatrix[nids][nsnps];
64
 
        char **gtMatrix = new (std::nothrow) char*[nids];
65
 
        for (unsigned int i=0;i<nids;i++) gtMatrix[i] = new (std::nothrow) char[nsnps];
66
 
 
67
 
        //Rprintf("1\n");
68
 
        std::string* Genotype;
69
 
        std::string sep="/";
70
 
        int nbytes;
71
 
 
72
 
        //Rprintf("nsnps=%d\n",nsnps);
73
 
        //Rprintf("nids=%d\n",nids);
74
 
 
75
 
        if ((nids % 4) == 0) nbytes = nidsTotal/4; else nbytes = ceil(1.*nidsTotal/4.);
76
 
 
77
 
        if (plink) sep=" ";
78
 
 
79
 
        if (append)
80
 
                fileWoA.open(filename.c_str(),std::fstream::app);
81
 
        else
82
 
                fileWoA.open(filename.c_str(),std::fstream::trunc);
83
 
 
84
 
        //Rprintf("A\n");
85
 
        for (unsigned int csnp=0;csnp<nsnps;csnp++) {
86
 
                // collect SNP data
87
 
                get_snps_many(snpdata+nbytes*csnp, &nidsTotal, &ieq1, gtint);
88
 
                for (int iii=from-1;iii<to;iii++) {
89
 
                        //Rprintf(" %d",gtint[iii]);
90
 
                        gtMatrix[iii-from+1][csnp]=gtint[iii];
91
 
                }
92
 
                //Rprintf("\n");
93
 
        }
94
 
        //Rprintf("B\n");
95
 
        for (int i=0;i<nids;i++) {
96
 
                fileWoA << i+from << " " << ids[i] << " 0 0 " << sex[i];
97
 
 
98
 
                for (int j=0;j<ntraits;j++) fileWoA << " " << 0;
99
 
                // unwrap genotypes
100
 
                for (unsigned int csnp=0;csnp<nsnps;csnp++) {
101
 
                        Genotype = getGenotype(coding[csnp],sep);
102
 
                        // figure out the coding
103
 
                        fileWoA << " " << Genotype[gtMatrix[i][csnp]];
104
 
                        //fileWoA << " x" << Letter0 << Letter1 << Genotype[0] << Genotype[1] << Genotype[2] << Genotype[3];
105
 
                }
106
 
                // end unwrap
107
 
                fileWoA << "\n";
108
 
        }
109
 
        //Rprintf("C\n");
110
 
        fileWoA.close();
111
 
 
112
 
        //Rprintf("oooo!" );
113
 
        //for (int i=0;i<10;i++) Rprintf("%d ",sex[i]);
114
 
        //Rprintf("oooo!\n" );
115
 
 
116
 
        sex.clear();
117
 
        delete [] gtMatrix;
118
 
        delete [] Genotype;
119
 
        delete [] gtint;
120
 
 
121
 
        return R_NilValue;
122
 
}
 
8
std::string* getGenotype(std::string coding, std::string sep)
 
9
{
 
10
    std::string* Genotype = new (std::nothrow) std::string [4];
 
11
    std::string Letter0 = coding.substr(0,1);
 
12
    std::string Letter1 = coding.substr(1,1);
 
13
    Genotype[0] = "0"+sep+"0";
 
14
    Genotype[1] = Letter0+sep+Letter0;
 
15
    Genotype[2] = Letter0+sep+Letter1;
 
16
    Genotype[3] = Letter1+sep+Letter1;
 
17
    return Genotype;
 
18
}
 
19
 
 
20
SEXP export_plink(SEXP Ids, SEXP Snpdata, SEXP Nsnps, SEXP NidsTotal,
 
21
                  SEXP Coding, SEXP From, SEXP To, SEXP Male, SEXP Traits,
 
22
                  SEXP Pedfilename, SEXP Plink, SEXP Append)
 
23
{
 
24
    int from = INTEGER(From)[0];
 
25
    int to = INTEGER(To)[0];
 
26
    
 
27
                
 
28
                if(from <1 || from > to) {error("The function SEXP export_plink(SEXP Ids, SEXP Snpdata, SEXP Nsnps, SEXP NidsTotal,... reports: the variable FROM should be >=1 and less then the variable TO.");} //Maksim 
 
29
 
 
30
                
 
31
                std::vector<unsigned short int> sex;
 
32
    sex.clear();
 
33
    unsigned short int sx;
 
34
    for(int i=(from - 1); i<to; i++) {
 
35
        sx = INTEGER(Male)[i];
 
36
        if (sx==0) sx=2;
 
37
        //Rprintf("%d %d\n",i,sx);
 
38
        sex.push_back(sx);
 
39
    }
 
40
    std::vector<std::string> ids;
 
41
    for(unsigned int i=0; i<((unsigned int) length(Ids)); i++)
 
42
        ids.push_back(CHAR(STRING_ELT(Ids,i)));
 
43
 
 
44
    std::vector<std::string> coding;
 
45
    for(unsigned int i=0; i<((unsigned int) length(Coding)); i++)
 
46
        coding.push_back(CHAR(STRING_ELT(Coding,i)));
 
47
 
 
48
    //Rprintf("0\n");
 
49
    unsigned int nsnps = INTEGER(Nsnps)[0];
 
50
    int nids = to - from + 1;
 
51
    int nidsTotal = INTEGER(NidsTotal)[0];
 
52
    int ntraits = INTEGER(Traits)[0];
 
53
    bool append = LOGICAL(Append)[0];
 
54
    bool plink = LOGICAL(Plink)[0];
 
55
    std::string filename = CHAR(STRING_ELT(Pedfilename,0));
 
56
    std::ofstream fileWoA;
 
57
    int ieq1 = 1;
 
58
    char * snpdata = (char *) RAW(Snpdata);
 
59
 
 
60
    //  int gtint[nidsTotal];
 
61
    int *gtint = new (std::nothrow) int[nidsTotal];
 
62
 
 
63
    //Rprintf("nsnps=%d\n",nsnps);
 
64
    //Rprintf("nids=%d\n",nids);
 
65
    //Rprintf("to=%d\n", to);
 
66
    //Rprintf("from=%d\n", from);
 
67
 
 
68
    //char gtMatrix[nids][nsnps];
 
69
    char **gtMatrix = new (std::nothrow) char*[nids];
 
70
    for (int i=0; i<nids; i++) {
 
71
        gtMatrix[i] = new (std::nothrow) char[nsnps];
 
72
    }
 
73
 
 
74
    //Rprintf("1\n");
 
75
    std::string* Genotype;
 
76
    std::string sep="/";
 
77
    int nbytes;
 
78
 
 
79
    //Rprintf("nsnps=%d\n",nsnps);
 
80
    //Rprintf("nids=%d\n",nids);
 
81
 
 
82
    if ((nids % 4) == 0) {
 
83
        nbytes = nidsTotal/4;
 
84
    }
 
85
    else {
 
86
        nbytes = ceil(1.*nidsTotal/4.);
 
87
    }
 
88
 
 
89
    if (plink) sep=" ";
 
90
 
 
91
    if (append)
 
92
        fileWoA.open(filename.c_str(),std::fstream::app);
 
93
    else
 
94
        fileWoA.open(filename.c_str(),std::fstream::trunc);
 
95
 
 
96
    //Rprintf("A\n");
 
97
    for (unsigned int csnp=0; csnp<nsnps; csnp++) {
 
98
        // collect SNP data
 
99
        get_snps_many(snpdata+nbytes*csnp, &nidsTotal, &ieq1, gtint);
 
100
        for (int iii=from-1; iii<to; iii++) {
 
101
            //Rprintf(" %d",gtint[iii]);
 
102
            gtMatrix[iii-from+1][csnp] = gtint[iii];
 
103
        }
 
104
        //Rprintf("\n");
 
105
    }
 
106
 
 
107
    //Rprintf("B\n");
 
108
    for (int i=0; i<nids; i++) {
 
109
        fileWoA << i+from << " " << ids[i] << " 0 0 " << sex[i];
 
110
 
 
111
        for (int j=0; j<ntraits; j++) fileWoA << " " << 0;
 
112
        // unwrap genotypes
 
113
        for (unsigned int csnp=0; csnp<nsnps; csnp++) {
 
114
            Genotype = getGenotype(coding[csnp], sep);
 
115
            // figure out the coding
 
116
            fileWoA << " " << Genotype[gtMatrix[i][csnp]];
 
117
            delete [] Genotype;
 
118
        }
 
119
        // end unwrap
 
120
        fileWoA << "\n";
 
121
    }
 
122
    //Rprintf("C\n");
 
123
    fileWoA.close();
 
124
 
 
125
    //Rprintf("oooo!" );
 
126
    //for (int i=0; i<10; i++) Rprintf("%d ",sex[i]);
 
127
    //Rprintf("oooo!\n" );
 
128
 
 
129
    sex.clear();
 
130
 
 
131
    for(int i=0; i<nids; i++) {
 
132
        delete [] gtMatrix[i];
 
133
    }
 
134
    delete [] gtMatrix;
 
135
 
 
136
    delete [] gtint;
 
137
 
 
138
    return R_NilValue;
 
139
}
 
140
 
123
141
 
124
142
SEXP export_plink_tped(SEXP Snpnames, SEXP Chromosomes, SEXP Map,
125
 
                SEXP Snpdata, SEXP Nsnps, SEXP Nids, SEXP Coding, SEXP Pedfilename,
126
 
                SEXP ExportNumeric)
 
143
                       SEXP Snpdata, SEXP Nsnps, SEXP Nids, SEXP Coding,
 
144
                       SEXP Pedfilename, SEXP ExportNumeric)
127
145
{
128
 
        std::vector<std::string> snpName;
129
 
        for(unsigned int i=0;i<((unsigned int) length(Snpnames));i++)
130
 
                snpName.push_back(CHAR(STRING_ELT(Snpnames,i)));
131
 
 
132
 
        std::vector<std::string> coding;
133
 
        for(unsigned int i=0;i<((unsigned int) length(Coding));i++)
134
 
                coding.push_back(CHAR(STRING_ELT(Coding,i)));
135
 
 
136
 
        std::vector<std::string> chromosome;
137
 
        for(unsigned int i=0;i<((unsigned int) length(Chromosomes));i++)
138
 
                chromosome.push_back(CHAR(STRING_ELT(Chromosomes,i)));
139
 
 
140
 
        std::vector<double> position;
141
 
        for(unsigned int i=0;i<((unsigned int) length(Map));i++)
142
 
                position.push_back(REAL(Map)[i]);
143
 
 
144
 
        //Rprintf("0\n");
145
 
        unsigned int nsnps = INTEGER(Nsnps)[0];
146
 
        int nids = INTEGER(Nids)[0];
147
 
        bool exportNumeric = LOGICAL(ExportNumeric)[0];
148
 
        std::string filename = CHAR(STRING_ELT(Pedfilename,0));
149
 
        std::ofstream fileWoA;
150
 
        int ieq1 = 1;
151
 
        char * snpdata = (char *) RAW(Snpdata);
152
 
 
153
 
        //      int gtint[nids];
154
 
        int *gtint = new (std::nothrow) int[nids];
155
 
 
156
 
        //Rprintf("nsnps=%d\n",nsnps);
157
 
        //Rprintf("nids=%d\n",nids);
158
 
 
159
 
        //Rprintf("1\n");
160
 
        std::string* Genotype;
161
 
        std::string sep=" ";
162
 
        int nbytes;
163
 
 
164
 
        //Rprintf("nsnps=%d\n",nsnps);
165
 
        //Rprintf("nids=%d\n",nids);
166
 
 
167
 
        if ((nids % 4) == 0) nbytes = nids/4; else nbytes = ceil(1.*nids/4.);
168
 
 
169
 
        fileWoA.open(filename.c_str(),std::fstream::trunc);
170
 
 
171
 
        //Rprintf("A\n");
172
 
        for (unsigned int csnp=0;csnp<nsnps;csnp++) {
173
 
                // collect SNP data
174
 
                get_snps_many(snpdata+nbytes*csnp, &nids, &ieq1, gtint);
175
 
                Genotype = getGenotype(coding[csnp],sep);
176
 
                fileWoA << chromosome[csnp] << " " << snpName[csnp] << " 0 " << (unsigned long int) position[csnp];
177
 
                if (!exportNumeric) {
178
 
                        for (int i=0;i<nids;i++) {
179
 
                                fileWoA << " " << Genotype[gtint[i]];
180
 
                        }
181
 
                } else {
182
 
                        for (int i=0;i<nids;i++) {
183
 
                                if (gtint[i]==0)
184
 
                                        fileWoA << " NA";
185
 
                                else
186
 
                                        fileWoA << " " << (gtint[i]-1);
187
 
                        }
188
 
                }
189
 
                fileWoA << "\n";
190
 
                //Rprintf("\n");
191
 
        }
192
 
        //Rprintf("B\n");
193
 
        /**
194
 
        for (int i=0;i<nids;i++) {
195
 
                fileWoA << i+from << " " << ids[i] << " 0 0 " << sex[i];
196
 
                for (int j=0;j<ntraits;j++) fileWoA << " " << 0;
197
 
                // unwrap genotypes
198
 
                for (unsigned int csnp=0;csnp<nsnps;csnp++) {
199
 
                        Genotype = getGenotype(coding[csnp],sep);
200
 
                        // figure out the coding
201
 
                        fileWoA << " " << Genotype[gtMatrix[i][csnp]];
202
 
                        //fileWoA << " x" << Letter0 << Letter1 << Genotype[0] << Genotype[1] << Genotype[2] << Genotype[3];
203
 
                }
204
 
                // end unwrap
205
 
                fileWoA << "\n";
206
 
        }
207
 
         **/
208
 
        //Rprintf("C\n");
209
 
        fileWoA.close();
210
 
 
211
 
        //Rprintf("oooo!" );
212
 
        //for (int i=0;i<10;i++) Rprintf("%d ",sex[i]);
213
 
        //Rprintf("oooo!\n" );
214
 
 
215
 
        delete [] Genotype;
216
 
        delete [] gtint;
217
 
 
218
 
        return R_NilValue;
 
146
    std::vector<std::string> snpName;
 
147
    for(unsigned int i=0; i<((unsigned int) length(Snpnames)); i++)
 
148
        snpName.push_back(CHAR(STRING_ELT(Snpnames,i)));
 
149
 
 
150
    std::vector<std::string> coding;
 
151
    for(unsigned int i=0; i<((unsigned int) length(Coding)); i++)
 
152
        coding.push_back(CHAR(STRING_ELT(Coding,i)));
 
153
 
 
154
    std::vector<std::string> chromosome;
 
155
    for(unsigned int i=0; i<((unsigned int) length(Chromosomes)); i++)
 
156
        chromosome.push_back(CHAR(STRING_ELT(Chromosomes,i)));
 
157
 
 
158
    std::vector<double> position;
 
159
    for(unsigned int i=0; i<((unsigned int) length(Map)); i++)
 
160
        position.push_back(REAL(Map)[i]);
 
161
 
 
162
    //Rprintf("0\n");
 
163
    unsigned int nsnps = INTEGER(Nsnps)[0];
 
164
    int nids = INTEGER(Nids)[0];
 
165
    bool exportNumeric = LOGICAL(ExportNumeric)[0];
 
166
    std::string filename = CHAR(STRING_ELT(Pedfilename,0));
 
167
    std::ofstream fileWoA;
 
168
    int ieq1 = 1;
 
169
    char * snpdata = (char *) RAW(Snpdata);
 
170
 
 
171
    //  int gtint[nids];
 
172
    int *gtint = new (std::nothrow) int[nids];
 
173
 
 
174
    //Rprintf("nsnps=%d\n",nsnps);
 
175
    //Rprintf("nids=%d\n",nids);
 
176
 
 
177
    //Rprintf("1\n");
 
178
    std::string* Genotype;
 
179
    std::string sep=" ";
 
180
    int nbytes;
 
181
 
 
182
    //Rprintf("nsnps=%d\n",nsnps);
 
183
    //Rprintf("nids=%d\n",nids);
 
184
 
 
185
    if ((nids % 4) == 0) {
 
186
        nbytes = nids/4;
 
187
    }
 
188
    else {
 
189
        nbytes = ceil(1.*nids/4.);
 
190
    }
 
191
 
 
192
    fileWoA.open(filename.c_str(), std::fstream::trunc);
 
193
 
 
194
    //Rprintf("A\n");
 
195
    for (unsigned int csnp=0; csnp<nsnps; csnp++) {
 
196
        // collect SNP data
 
197
        get_snps_many(snpdata+nbytes*csnp, &nids, &ieq1, gtint);
 
198
        Genotype = getGenotype(coding[csnp], sep);
 
199
        fileWoA << chromosome[csnp] << " " << snpName[csnp]
 
200
                << " 0 " << (unsigned long int) position[csnp];
 
201
 
 
202
        if (!exportNumeric) {
 
203
            for (int i=0; i<nids; i++) {
 
204
                fileWoA << " " << Genotype[gtint[i]];
 
205
            }
 
206
        } else {
 
207
            for (int i=0; i<nids; i++) {
 
208
                if (gtint[i]==0)
 
209
                    fileWoA << " NA";
 
210
                else
 
211
                    fileWoA << " " << (gtint[i]-1);
 
212
            }
 
213
        }
 
214
        fileWoA << "\n";
 
215
    delete [] Genotype;
 
216
        //Rprintf("\n");
 
217
    }
 
218
    //Rprintf("B\n");
 
219
    /**
 
220
       for (int i=0; i<nids; i++) {
 
221
       fileWoA << i+from << " " << ids[i] << " 0 0 " << sex[i];
 
222
       for (int j=0; j<ntraits; j++) fileWoA << " " << 0;
 
223
       // unwrap genotypes
 
224
       for (unsigned int csnp=0; csnp<nsnps; csnp++) {
 
225
       Genotype = getGenotype(coding[csnp],sep);
 
226
       // figure out the coding
 
227
       fileWoA << " " << Genotype[gtMatrix[i][csnp]];
 
228
       //fileWoA << " x" << Letter0 << Letter1 << Genotype[0] << Genotype[1] << Genotype[2] << Genotype[3];
 
229
       }
 
230
       // end unwrap
 
231
       fileWoA << "\n";
 
232
       }
 
233
    **/
 
234
    //Rprintf("C\n");
 
235
    fileWoA.close();
 
236
 
 
237
    //Rprintf("oooo!" );
 
238
    //for (int i=0; i<10; i++) Rprintf("%d ",sex[i]);
 
239
    //Rprintf("oooo!\n" );
 
240
 
 
241
    delete [] gtint;
 
242
 
 
243
    return R_NilValue;
219
244
}
220
245
 
221
246
#ifdef __cplusplus