ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/tophat_cpp/long_spanning_reads.cpp
(Generate patch)
# Line 1179 | Line 1179
1179        vector<FILE*>& possible_juncs_files,
1180        vector<FILE*>& possible_insertions_files,
1181        vector<FILE*>& possible_deletions_files,
1182 <      vector<FZPipe>& spliced_seg_files,
1183 <      vector<FZPipe>& seg_files,
1182 >      vector<FZPipe>& spliced_segmap_files, //.sam.z files
1183 >      vector<string>& segmap_fnames, //.bam files
1184        FZPipe& reads_file)
1185   {
1186 <  if (seg_files.size() == 0)
1186 >  if (segmap_fnames.size() == 0)
1187    {
1188      fprintf(stderr, "No hits to process, exiting\n");
1189      exit(0);
# Line 1205 | Line 1205
1205    vector<HitStream> spliced_hits;
1206  
1207    vector<HitFactory*> factories;
1208 <  for (size_t i = 0; i < seg_files.size(); ++i)
1208 >  for (size_t i = 0; i < segmap_fnames.size(); ++i)
1209    {
1210 <    HitFactory* fac = new BowtieHitFactory(it, rt);
1210 >    HitFactory* fac = new BAMHitFactory(it, rt);
1211      factories.push_back(fac);
1212 <    HitStream hs(seg_files[i].file, fac, false, false, false, need_seq, need_qual);
1212 >    HitStream hs(segmap_fnames[i], fac, false, false, false, need_seq, need_qual);
1213      contig_hits.push_back(hs);
1214    }
1215  
1216    fprintf(stderr, "Loading spliced hits...");
1217  
1218    //vector<vector<pair<uint32_t, HitsForRead> > > spliced_hits;
1219 <  for (size_t i = 0; i < spliced_seg_files.size(); ++i)
1219 >  for (size_t i = 0; i < spliced_segmap_files.size(); ++i)
1220    {
1221      int anchor_length = 0;
1222  
# Line 1224 | Line 1224
1224      // if (i == 0 || i == spliced_seg_files.size() - 1)
1225      // anchor_length = min_anchor_len;
1226  
1227 <    HitFactory* fac = new SplicedBowtieHitFactory(it,
1227 >    HitFactory* fac = new SplicedSAMHitFactory(it,
1228                      rt,
1229                      anchor_length);
1230      factories.push_back(fac);
1231  
1232 <    HitStream hs(spliced_seg_files[i].file, fac, true, false, false, need_seq, need_qual);
1232 >    HitStream hs(spliced_segmap_files[i].file, fac, true, false, false, need_seq, need_qual);
1233      spliced_hits.push_back(hs);
1234  
1235    }
# Line 1327 | Line 1327
1327    join_segment_hits(bam_writer, possible_juncs, possible_insertions, it, rt, reads_file.file, contig_hits, spliced_hits);
1328    //join_segment_hits(possible_juncs, possible_insertions, it, rt, reads_file.file, contig_hits, spliced_hits);
1329  
1330 <  for (size_t i = 0; i < seg_files.size(); ++i)
1331 <      seg_files[i].close();
1332 <    for (size_t i = 0; i < spliced_seg_files.size(); ++i)
1333 <        spliced_seg_files[i].close();
1330 >  /*for (size_t i = 0; i < segmap_fnames.size(); ++i)
1331 >      segmap_fnames[i].close();
1332 >  */
1333 >  for (size_t i = 0; i < spliced_segmap_files.size(); ++i)
1334 >        spliced_segmap_files[i].close();
1335      reads_file.close();
1336  
1337    for (size_t fac = 0; fac < factories.size(); ++fac)
# Line 1395 | Line 1396
1396        }
1397  
1398  
1399 <  string segment_file_list = argv[optind++];
1400 <  string spliced_segment_file_list;
1399 >  string segmap_file_list = argv[optind++];
1400 >  string spliced_segmap_flist;
1401    if(optind < argc)
1402      {
1403 <      spliced_segment_file_list = argv[optind++];
1403 >      spliced_segmap_flist = argv[optind++];
1404      }
1405  
1406    ifstream ref_stream(ref_file_name.c_str(), ifstream::in);
# Line 1409 | Line 1410
1410    checkSamHeader();
1411  
1412    //FILE* reads_file = fopen(reads_file_name.c_str(), "r");
1413 <  vector<string> segment_file_names;
1414 <  tokenize(segment_file_list, ",",segment_file_names);
1413 >  vector<string> segmap_file_names;
1414 >  tokenize(segmap_file_list, ",",segmap_file_names);
1415  
1416 <  string unzcmd=getUnpackCmd(reads_file_name, spliced_segment_file_list.empty() &&
1417 <                 segment_file_names.size()<4);
1416 >  string unzcmd=getUnpackCmd(reads_file_name, spliced_segmap_flist.empty() &&
1417 >                 segmap_file_names.size()<4);
1418    FZPipe reads_file(reads_file_name, unzcmd);
1419    if (reads_file.file==NULL)
1420       err_die("Error: cannot open %s for reading\n",
# Line 1475 | Line 1476
1476          insertions_files.push_back(insertions_file);
1477      }
1478  
1479 <    //vector<FILE*> segment_files;
1480 <    vector<FZPipe> segment_files;
1481 <    for (size_t i = 0; i < segment_file_names.size(); ++i)
1482 <    {
1483 <      fprintf(stderr, "Opening %s for reading\n",
1484 <        segment_file_names[i].c_str());
1485 <      //FILE* seg_file = fopen(segment_file_names[i].c_str(), "r");
1486 <      FZPipe seg_file(segment_file_names[i], unzcmd);
1479 >    //vector<FZPipe> segmap_files;
1480 >    /*
1481 >    vector<string> segmap_files; //BAM files
1482 >    for (size_t i = 0; i < segmap_file_names.size(); ++i)
1483 >    {
1484 >        fprintf(stderr, "Opening %s for reading\n",
1485 >        segmap_file_names[i].c_str());
1486 >      FZPipe seg_file(segmap_file_names[i], unzcmd);
1487        if (seg_file.file == NULL)
1488          {
1489          err_die("Error: cannot open %s for reading\n",
1490 <              segment_file_names[i].c_str());
1490 >              segmap_file_names[i].c_str());
1491          }
1492 <      segment_files.push_back(seg_file);
1492 <    }
1492 >      segmap_files.push_back(seg_file);
1493  
1494 <  vector<string> spliced_segment_file_names;
1495 <  //vector<FILE*> spliced_segment_files;
1496 <  vector<FZPipe> spliced_segment_files;
1497 <  tokenize(spliced_segment_file_list, ",",spliced_segment_file_names);
1498 <  for (size_t i = 0; i < spliced_segment_file_names.size(); ++i)
1499 <    {
1500 <      fprintf(stderr, "Opening %s for reading\n",
1501 <        spliced_segment_file_names[i].c_str());
1502 <      //FILE* spliced_seg_file = fopen(spliced_segment_file_names[i].c_str(), "r");
1503 <      FZPipe spliced_seg_file(spliced_segment_file_names[i], unzcmd);
1494 >    }*/
1495 >
1496 >  vector<string> spliced_segmap_file_names;
1497 >  vector<FZPipe> spliced_segmap_files;
1498 >  //these are headerless .sam.z files
1499 >  tokenize(spliced_segmap_flist, ",",spliced_segmap_file_names);
1500 >  for (size_t i = 0; i < spliced_segmap_file_names.size(); ++i)
1501 >    {
1502 >      if (i==0) {
1503 >        unzcmd=getUnpackCmd(spliced_segmap_file_names[i]);
1504 >        }
1505 >          fprintf(stderr, "Opening %s for reading\n",
1506 >        spliced_segmap_file_names[i].c_str());
1507 >      FZPipe spliced_seg_file(spliced_segmap_file_names[i], unzcmd);
1508        if (spliced_seg_file.file == NULL)
1509          {
1510          err_die("Error: cannot open %s for reading\n",
1511 <             spliced_segment_file_names[i].c_str());
1511 >             spliced_segmap_file_names[i].c_str());
1512          }
1513 <        spliced_segment_files.push_back(spliced_seg_file);
1513 >        spliced_segmap_files.push_back(spliced_seg_file);
1514      }
1515  
1516 <  if (spliced_segment_files.size())
1516 >  if (spliced_segmap_files.size())
1517      {
1518 <      if (spliced_segment_files.size() != segment_files.size())
1518 >      if (spliced_segmap_files.size() != segmap_file_names.size())
1519          {
1520 <        err_die("Error: each segment file must have a corresponding spliced segment file\n");
1520 >        err_die("Error: each segment mapping file must have a corresponding spliced segment mapping file\n");
1521          }
1522      }
1523    GBamWriter bam_writer("-", sam_header.c_str());
1524 <  driver(bam_writer, ref_stream, juncs_files, insertions_files, deletions_files, spliced_segment_files, segment_files, reads_file);
1521 <  //driver(ref_stream, juncs_files, insertions_files, deletions_files, spliced_segment_files, segment_files, reads_file);
1524 >  driver(bam_writer, ref_stream, juncs_files, insertions_files, deletions_files, spliced_segmap_files, segmap_file_names, reads_file);
1525    return 0;
1526   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines