Commit c427ed87 authored by Li's avatar Li
Browse files

Output mapq in BED format. Update readme

parent 90246ae4
Loading
Loading
Loading
Loading
+8 −1
Original line number Diff line number Diff line
@@ -113,6 +113,13 @@ threshold to make a correction by setting **--bc-probability-threshold**
corrections. For scATAC-seq data with multiple read and barcode files, you can
use "," to concatenate multiple input files as the example [above](#general). 

The BED format (fragment file) for bulk and single-cell are different except for the first
three columns. For bulk data, the columns are "chrom chrom_start chrom_end N mapq strand".
For single-cell data, the columns are "chrom chrom_start chrom_end barcode duplicate_count" 
as the definition of the fragment file in 
[CellRanger](https://support.10xgenomics.com/single-cell-atac/software/pipelines/latest/output/fragments). 
Note that chrom_end is open-end.

#### <a name="map-hic"></a>Map Hi-C short reads

```sh
+5 −5
Original line number Diff line number Diff line
@@ -1033,7 +1033,7 @@ class BEDOutputTools : public OutputTools<MappingRecord> {
    std::string strand = mapping.IsPositive() ? "+" : "-";
    const char *reference_sequence_name = reference.GetSequenceNameAt(rid);
    uint32_t mapping_end_position = mapping.GetEndPosition();
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t1000\t" + strand + "\n");
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t" + std::to_string(mapping.mapq) + "\t" + strand + "\n");
  }
};

@@ -1053,7 +1053,7 @@ class BEDPEOutputTools : public OutputTools<MappingRecord> {
    std::string strand = mapping.IsPositive() ? "+" : "-";
    const char *reference_sequence_name = reference.GetSequenceNameAt(rid);
    uint32_t mapping_end_position = mapping.GetEndPosition();
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t1000\t" + strand + "\n");
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t" + std::to_string(mapping.mapq) + "\t" + strand + "\n");
  }
};

@@ -1075,7 +1075,7 @@ class TagAlignOutputTools : public OutputTools<MappingRecord> {
    std::string strand = mapping.IsPositive() ? "+" : "-";
    const char *reference_sequence_name = reference.GetSequenceNameAt(rid);
    uint32_t mapping_end_position = mapping.GetEndPosition();
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t1000\t" + strand + "\n");
    this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping_end_position) + "\tN\t" + std::to_string(mapping.mapq) + "\t" + strand + "\n");
  }
};

@@ -1090,9 +1090,9 @@ class PairedTagAlignOutputTools : public OutputTools<MappingRecord> {
    uint32_t negative_read_start = negative_read_end - mapping.negative_alignment_length;
    const char *reference_sequence_name = reference.GetSequenceNameAt(rid);
    if (positive_strand) {
      this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.fragment_start_position) + "\t" + std::to_string(positive_read_end) + "\tN\t1000\t+\n" + std::string(reference_sequence_name) + "\t" + std::to_string(negative_read_start) + "\t" + std::to_string(negative_read_end) + "\tN\t1000\t-\n");
      this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(mapping.fragment_start_position) + "\t" + std::to_string(positive_read_end) + "\tN\t" + std::to_string(mapping.mapq) + "\t+\n" + std::string(reference_sequence_name) + "\t" + std::to_string(negative_read_start) + "\t" + std::to_string(negative_read_end) + "\tN\t" + std::to_string(mapping.mapq) + "\t-\n");
    } else {
      this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(negative_read_start) + "\t" + std::to_string(negative_read_end) + "\tN\t1000\t-\n" + std::string(reference_sequence_name) + "\t" + std::to_string(mapping.fragment_start_position) + "\t" + std::to_string(positive_read_end) + "\tN\t1000\t+\n");
      this->AppendMappingOutput(std::string(reference_sequence_name) + "\t" + std::to_string(negative_read_start) + "\t" + std::to_string(negative_read_end) + "\tN\t" + std::to_string(mapping.mapq) + "\t-\n" + std::string(reference_sequence_name) + "\t" + std::to_string(mapping.fragment_start_position) + "\t" + std::to_string(positive_read_end) + "\tN\t" + std::to_string(mapping.mapq) + "\t+\n");
    }
  }
};