« back to all changes in this revision
Viewing changes to cluster.cpp
- browse files at revision 7
- compare with another revision
- download diff
- download tarball
- view history from revision 7
- Committer: Bazaar Package Importer
- Author(s): Andreas Tille
- Date: 2009-10-23 13:35:02 UTC
- mfrom: (1.1.2 upstream)
- Revision ID: james.westby@ubuntu.com-20091023133502-rmnkrfi314mevnvt
Tags: 1.07-1
* New upstream version
* Removed debian/patches/20_plink-1.06-gcc4.4.patch because it
is applied upstream
* Standards-Version: 3.8.3 (no changes needed)
* Debhelper 7
* Build-Depends: zlib1g-dev
* Removed debian/patches/20_plink-1.06-gcc4.4.patch because it
is applied upstream
* Standards-Version: 3.8.3 (no changes needed)
* Debhelper 7
* Build-Depends: zlib1g-dev
- files added:
- annot.cpp
- files removed:
- debian/patches/20_plink-1.06-gcc4.4.patch
- files modified:
- Makefile
added
removed
cluster.cpp
2
2
3
3
//////////////////////////////////////////////////////////////////
4
4
// //
5
// PLINK (c) 2005-2008 Shaun Purcell //
5
// PLINK (c) 2005-2009 Shaun Purcell //
6
6
// //
7
7
// This file is distributed under the GNU General Public //
8
8
// License, Version 2. Please see the file COPYING for more //
411
411
else // ... read IBS information from .genome file
412
412
{
413
413
414
415
414
416
checkFileExists(par::ibd_file);
415
417
416
418
if ( par::ibd_read_minimal )
421
423
+par::ibd_file+" ] \n");
422
424
423
425
424
ifstream INC;
425
INC.open(par::ibd_file.c_str());
426
if ( compressed( par::ibd_file ) )
427
par::compress_genome = true;
428
429
ZInput ZINC( par::ibd_file , par::compress_genome );
430
426
431
427
432
map<string,int> mperson;
428
433
for (int i=0; i<n; i++)
440
445
// read in list of people here
441
446
while ( 1 )
442
447
{
443
vector<string> ids = tokenizeLine(INC);
448
vector<string> ids = ZINC.tokenizeLine();
444
449
if ( ids.size() != 2 )
445
450
{
446
451
string emsg = "Problem with line in [ " + par::ibd_file + " ]\n";
472
477
473
478
// Just in case we have a malformed file
474
479
475
if ( INC.eof() )
480
if ( ZINC.endOfFile() )
476
481
error("Problem with premature stop in file [ " + par::ibd_file + " ]\n");
477
482
478
483
}
479
484
485
480
486
//////////////////////////////////////////////////////
481
487
// Now read the actual IBS/PPC values for these peeps
482
488
491
497
{
492
498
double mydst, pv, ibd;
493
499
494
vector<string> val = tokenizeLine(INC);
500
vector<string> val = ZINC.tokenizeLine();
495
501
496
if ( INC.eof() )
502
if ( ZINC.endOfFile() )
497
503
{
498
504
// Check that p1,p2 counts are as should be...
499
505
break;
606
612
int col_length = 0;
607
613
double mydst;
608
614
609
vector<string> tokens = tokenizeLine(INC);
615
vector<string> tokens = ZINC.tokenizeLine();
610
616
col_length = tokens.size();
611
617
612
618
if ( tokens.size() < 4 ||
631
637
error("Could not find PPC or DST fields in .genome file");
632
638
633
639
// Read each pair at a time
634
while ( ! INC.eof() )
640
while ( ! ZINC.endOfFile() )
635
641
{
636
642
637
vector<string> tokens = tokenizeLine(INC);
643
vector<string> tokens = ZINC.tokenizeLine();
638
644
639
645
if ( tokens.size() == 0 )
640
646
continue;
654
660
string ipv = tokens[ppc_code];
655
661
string idst = tokens[dst_code];
656
662
663
// Skip any blank rows, or additional header rows
657
664
if (fid1=="") continue;
665
if (fid1=="FID1") continue;
658
666
659
667
// if ( ! ( from_string<double>( ibs0 , i0 , std::dec) &&
660
668
// from_string<double>( ibs1 , i1 , std::dec) &&
722
730
}
723
731
724
732
725
INC.close();
733
ZINC.close();
726
734
727
735
728
736
/////////////////////////////////////////////
1568
1576
1569
1577
}
1570
1578
1579
1580
void Plink::groupGenome()
1581
{
1582
1583
// Read from a (non-verbose) genome file
1584
1585
checkFileExists(par::ibd_file);
1586
1587
if ( par::ibd_read_minimal )
1588
printLOG("Reading IBS estimates (minimal format) from [ "
1589
+par::ibd_file+" ] \n");
1590
else
1591
printLOG("Reading genome-wide IBS estimates from [ "
1592
+par::ibd_file+" ] \n");
1593
1594
ifstream INC;
1595
INC.open(par::ibd_file.c_str());
1596
1597
1598
map<string,int> mperson;
1599
for (int i=0; i<n; i++)
1600
mperson.insert(make_pair( sample[i]->fid+"_"+sample[i]->iid , i ));
1601
1602
map<Individual*,int> mcode;
1603
for (int i=0; i<n; i++)
1604
mcode.insert(make_pair( sample[i] , i ));
1605
1606
vector<Individual*> peeps;
1607
1608
// We wish to read in an NxN matrix, and convert it to a KxK one
1609
1610
matrix_t dk( nk );
1611
sizeMatrix(dk,nk,0);
1612
for (int j=0; j<nk; j++)
1613
dk[j].resize(j,0);
1614
1615
table_t dkn;
1616
sizeTable(dkn,nk,0);
1617
for (int j=0; j<nk; j++)
1618
dkn[j].resize(j,0);
1619
1620
1621
// Read in .genome file in verbose mode
1622
1623
// We only want FID1,IID1,FID2,IID2 (always first four)
1624
// DST and PPC
1625
1626
// Get field codes from header
1627
1628
int dst_code = -1;
1629
int col_length = 0;
1630
double mydst;
1631
1632
vector<string> tokens = tokenizeLine(INC);
1633
col_length = tokens.size();
1634
1635
if ( tokens.size() < 4 ||
1636
tokens[0] != "FID1" ||
1637
tokens[1] != "IID1" ||
1638
tokens[2] != "FID2" ||
1639
tokens[3] != "IID2" )
1640
error("Problem with header row of .genome file");
1641
1642
1643
for ( int i = 4; i<tokens.size(); i++)
1644
{
1645
if ( tokens[i] == "DST" )
1646
dst_code = i;
1647
}
1648
1649
if ( dst_code == -1 )
1650
error("Could not find DST fields in .genome file");
1651
1652
// Read each pair at a time
1653
while ( ! INC.eof() )
1654
{
1655
1656
vector<string> tokens = tokenizeLine(INC);
1657
1658
if ( tokens.size() == 0 )
1659
continue;
1660
1661
if ( col_length != tokens.size() )
1662
{
1663
string strmsg = "";
1664
for (int i=0;i<tokens.size();i++)
1665
strmsg += tokens[i] + " ";
1666
error("Problem reading line in .genome file:\n"+strmsg+"\n");
1667
}
1668
1669
string fid1 = tokens[0];
1670
string iid1 = tokens[1];
1671
string fid2 = tokens[2];
1672
string iid2 = tokens[3];
1673
string idst = tokens[dst_code];
1674
1675
if (fid1=="") continue;
1676
1677
if ( ! from_string<double>( mydst , idst , std::dec) )
1678
mydst = 0;
1679
1680
map<string,int>::iterator person1 = mperson.find(fid1+"_"+iid1);
1681
map<string,int>::iterator person2 = mperson.find(fid2+"_"+iid2);
1682
1683
if ( person1 == mperson.end() ||
1684
person2 == mperson.end() ||
1685
person1 == person2 )
1686
continue;
1687
1688
int k1 = sample[ person1->second ]->sol;
1689
int k2 = sample[ person2->second ]->sol;
1690
1691
if ( k1 < 0 || k2 < 0 || k1 == k2 )
1692
continue;
1693
1694
if ( k2 > k1 )
1695
{
1696
int tmp = k2;
1697
k2 = k1;
1698
k1 = tmp;
1699
}
1700
1701
// Record IBS distance
1702
1703
dk[k1][k2] = mydst;
1704
++dkn[k1][k2];
1705
1706
} // Read next line in .genome
1707
1708
INC.close();
1709
1710
for (int i=0; i<nk; i++)
1711
for (int j=0; j<i; j++)
1712
{
1713
if ( dkn[i][j] > 0 )
1714
dk[i][j] /= (double)dkn[i][j];
1715
}
1716
1717
// Output a dummy .genome file
1718
1719
ofstream GOUT0;
1720
GOUT0.open( (par::output_file_name + ".plst").c_str(), ios::out);
1721
printLOG("Writing person include list to [ " + par::output_file_name + ".plst ]\n");
1722
1723
ofstream GOUT1;
1724
GOUT1.open( (par::output_file_name + ".clst").c_str(), ios::out);
1725
printLOG("Writing cluster list to [ " + par::output_file_name + ".clst ]\n");
1726
1727
1728
map<int,int> k2i;
1729
for (int i=0;i<n; i++)
1730
{
1731
int j = sample[i]->sol;
1732
if ( j < 0 )
1733
continue;
1734
if ( k2i.find(j) == k2i.end() )
1735
{
1736
k2i.insert(make_pair(j,i));
1737
GOUT0 << sample[i]->fid << " "
1738
<< sample[i]->iid << "\n";
1739
GOUT1 << sample[i]->fid << " "
1740
<< sample[i]->iid << " "
1741
<< kname[j] << "\n";
1742
}
1743
}
1744
GOUT0.close();
1745
GOUT1.close();
1746
1747
ofstream GOUT;
1748
GOUT.open( (par::output_file_name + ".genome").c_str(), ios::out);
1749
printLOG("Writing grouped .genome file to [ " + par::output_file_name + ".genome ]\n");
1750
GOUT << setw(par::pp_maxfid) << "FID1" << " "
1751
<< setw(par::pp_maxiid) << "IID1" << " "
1752
<< setw(par::pp_maxfid) << "FID2" << " "
1753
<< setw(par::pp_maxiid) << "IID2" << " "
1754
<< setw(8) << "DST" << " "
1755
<< setw(8) << "PPC" << "\n";
1756
1757
for (int i=0; i<nk; i++)
1758
for (int j=0; j<i; j++)
1759
{
1760
Individual * s1 = sample[k2i.find(i)->second];
1761
Individual * s2 = sample[k2i.find(j)->second];
1762
1763
GOUT << s1->fid << " " << s1->iid << " "
1764
<< s2->fid << " " << s2->iid << " ";
1765
1766
// cout << i << " " << j << " ";
1767
// cout << dk.size() << " " << dk[i].size() << "\n";
1768
1769
// cout << dkn.size() << " " << dkn[i].size() << "\n";
1770
// cout << dk[i][j] << " of " << dkn[i][j] << "\n";
1771
1772
GOUT << dk[i][j] << " 1\n";
1773
}
1774
1775
GOUT.close();
1776
1777
}
Loggerhead is a web-based interface for Breezy
Version: 2.0.1