Commit b75d1205 authored by Haowen Zhang's avatar Haowen Zhang
Browse files

Add CB tag for corrected barcodes in SAM

parent 2c14e3de

src/chromap.cc

+6 −6
Original line number Diff line number Diff line
@@ -1754,8 +1754,8 @@ void Chromap<MappingRecord>::ProcessBestMappingsForPairedEndReadOnOneDirection(D 
				  }

				//printf("%d %d\n", ref_start_position1, ref_end_position1);

				//printf("%d %d\n", ref_start_position2, ref_end_position2);

				  EmplaceBackMappingRecord(read_id, read1_name, 1, ref_start_position1, rid1, flag1, first_read_direction == kPositive ? 1 : 0, is_unique, mapq, NM1, n_cigar1, cigar1, MD_tag1, &((*mappings_on_diff_ref_seqs)[rid1])); 

				  EmplaceBackMappingRecord(read_id, read2_name, 1, ref_start_position2, rid2, flag2, second_read_direction == kPositive ? 1 : 0, is_unique, mapq, NM2, n_cigar2, cigar2, MD_tag2, &((*mappings_on_diff_ref_seqs)[rid2])); 

				  EmplaceBackMappingRecord(read_id, read1_name, barcode_key, 1, ref_start_position1, rid1, flag1, first_read_direction == kPositive ? 1 : 0, is_unique, mapq, NM1, n_cigar1, cigar1, MD_tag1, &((*mappings_on_diff_ref_seqs)[rid1])); 

				  EmplaceBackMappingRecord(read_id, read2_name, barcode_key, 1, ref_start_position2, rid2, flag2, second_read_direction == kPositive ? 1 : 0, is_unique, mapq, NM2, n_cigar2, cigar2, MD_tag2, &((*mappings_on_diff_ref_seqs)[rid2])); 

				} else if (output_mapping_in_pairs_) {

					int position1 = ref_start_position1;

					int position2 = ref_start_position2;
@@ -2123,12 +2123,12 @@ void Chromap<PAFMapping>::EmplaceBackMappingRecord(uint32_t read_id, const char 
}



template<typename MappingRecord>

void Chromap<MappingRecord>::EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<MappingRecord> *mappings_on_diff_ref_seqs) {

void Chromap<MappingRecord>::EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint32_t cell_barcode, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<MappingRecord> *mappings_on_diff_ref_seqs) {

}



template<>

void Chromap<SAMMapping>::EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<SAMMapping> *mappings_on_diff_ref_seqs) {

  mappings_on_diff_ref_seqs->emplace_back(SAMMapping{read_id, std::string(read_name), num_dups, position, rid, flag, direction, 0, is_unique, mapq, NM, n_cigar, cigar, MD_tag});

void Chromap<SAMMapping>::EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint32_t cell_barcode, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<SAMMapping> *mappings_on_diff_ref_seqs) {

  mappings_on_diff_ref_seqs->emplace_back(SAMMapping{read_id, std::string(read_name), cell_barcode, num_dups, position, rid, flag, direction, 0, is_unique, mapq, NM, n_cigar, cigar, MD_tag});

}



template <typename MappingRecord>
@@ -2425,7 +2425,7 @@ void Chromap<MappingRecord>::ProcessBestMappingsForSingleEndRead(Direction mappi 
		if (*num_best_mappings_reported >= 1) {

			flag |= BAM_FSECONDARY;

		}

		EmplaceBackMappingRecord(read_id, read_name, 1, ref_start_position, rid, flag, 0, is_unique, mapq, NM, n_cigar, cigar, MD_tag, &((*mappings_on_diff_ref_seqs)[rid])); 

		EmplaceBackMappingRecord(read_id, read_name, barcode_key, 1, ref_start_position, rid, flag, 0, is_unique, mapq, NM, n_cigar, cigar, MD_tag, &((*mappings_on_diff_ref_seqs)[rid])); 

	} else if (output_mapping_in_PAF_) {

        	EmplaceBackMappingRecord(read_id, read_name, read_length, barcode_key, ref_start_position, ref_end_position - ref_start_position + 1, mapq, direction, is_unique, 1, &((*mappings_on_diff_ref_seqs)[rid]));

	} else {

src/chromap.h

+1 −1
Original line number Diff line number Diff line
@@ -138,7 +138,7 @@ class Chromap { 
  void ProcessBestMappingsForSingleEndRead(Direction mapping_direction, uint8_t mapq, int num_candidates, uint32_t repetitive_seed_length, int min_num_errors, int num_best_mappings, int second_min_num_errors, int num_second_best_mappings, const SequenceBatch &read_batch, uint32_t read_index, const SequenceBatch &reference, const SequenceBatch &barcode_batch, const std::vector<int> &best_mapping_indices, const std::vector<std::pair<int, uint64_t> > &mappings, const std::vector<int> &split_sites, int *best_mapping_index, int *num_best_mappings_reported, std::vector<std::vector<MappingRecord> > *mappings_on_diff_ref_seqs);

  void EmplaceBackMappingRecord(uint32_t read_id, uint32_t barcode, uint32_t fragment_start_position, uint16_t fragment_length, uint8_t mapq, uint8_t direction, uint8_t is_unique, uint8_t num_dups, std::vector<MappingRecord> *mappings_on_diff_ref_seqs);

  void EmplaceBackMappingRecord(uint32_t read_id, const char* read_name, uint16_t read_length, uint32_t barcode, uint32_t fragment_start_position, uint16_t fragment_length, uint8_t mapq, uint8_t direction, uint8_t is_unique, uint8_t num_dups, std::vector<MappingRecord> *mappings_on_diff_ref_seqs);

  void EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<MappingRecord> *mappings_on_diff_ref_seqs);

  void EmplaceBackMappingRecord(uint32_t read_id, const char *read_name, uint32_t cell_barcode, uint8_t num_dups, int64_t position, int rid, int flag, uint8_t direction, uint8_t is_unique, uint8_t mapq, uint32_t NM, int n_cigar, uint32_t *cigar, std::string &MD_tag, std::vector<MappingRecord> *mappings_on_diff_ref_seqs);

  void ApplyTn5ShiftOnSingleEndMapping(uint32_t num_reference_sequences, std::vector<std::vector<MappingRecord> > *mappings);

  uint32_t LoadSingleEndReadsWithBarcodes(SequenceBatch *read_batch, SequenceBatch *barcode_batch);

src/output_tools.h

+7 −4
Original line number Diff line number Diff line
@@ -298,6 +298,7 @@ struct PairedPAFMapping { 
struct SAMMapping {

  uint32_t read_id;

  std::string read_name;

  uint32_t cell_barcode;

  //uint16_t read_length;

  //uint32_t fragment_start_position;

  //uint16_t fragment_length;
@@ -318,16 +319,16 @@ struct SAMMapping { 
  bool operator<(const SAMMapping& m) const {

    int read1_flag = flag & BAM_FREAD1;

    int m_read1_flag = m.flag & BAM_FREAD1;

    return std::tie(rid, pos, mapq, read_id, read1_flag) < std::tie(m.rid, m.pos, m.mapq, m.read_id, m_read1_flag);

    return std::tie(rid, pos, cell_barcode, mapq, read_id, read1_flag) < std::tie(m.rid, m.pos, m.cell_barcode, m.mapq, m.read_id, m_read1_flag);

  }

  bool operator==(const SAMMapping& m) const {

    return std::tie(pos, rid, is_rev) == std::tie(m.pos, m.rid, m.is_rev);

    return std::tie(pos, rid, cell_barcode, is_rev) == std::tie(m.pos, m.rid, m.cell_barcode, m.is_rev);

  }

  bool HasSamePosition(const SAMMapping& m) const {

    return std::tie(pos, rid, is_rev) == std::tie(m.pos, m.rid, m.is_rev);

  }

  uint32_t GetBarcode() const {

    return 0;

    return cell_barcode;

  }

  void Tn5Shift() {

    // We don't support Tn5 shift in SAM format because it has other fields that depend mapping position.
@@ -389,6 +390,7 @@ struct SAMMapping { 
    uint16_t read_name_length = read_name.length();

    num_written_bytes += fwrite(&read_name_length, sizeof(uint16_t), 1, temp_mapping_output_file);

    num_written_bytes += fwrite(read_name.data(), sizeof(char), read_name_length, temp_mapping_output_file);

    num_written_bytes += fwrite(&cell_barcode, sizeof(uint32_t), 1, temp_mapping_output_file);

    num_written_bytes += fwrite(&num_dups, sizeof(uint8_t), 1, temp_mapping_output_file);

    num_written_bytes += fwrite(&pos, sizeof(int64_t), 1, temp_mapping_output_file);

    num_written_bytes += fwrite(&rid, sizeof(int), 1, temp_mapping_output_file);
@@ -413,6 +415,7 @@ struct SAMMapping { 
    num_read_bytes += fread(&read_name_length, sizeof(uint16_t), 1, temp_mapping_output_file);

    read_name = std::string(read_name_length, '\0');

    num_read_bytes += fread(&(read_name[0]), sizeof(char), read_name_length, temp_mapping_output_file);

    num_read_bytes += fread(&cell_barcode, sizeof(uint32_t), 1, temp_mapping_output_file);

    num_read_bytes += fread(&num_dups, sizeof(uint8_t), 1, temp_mapping_output_file);

    num_read_bytes += fread(&pos, sizeof(int64_t), 1, temp_mapping_output_file);

    num_read_bytes += fread(&rid, sizeof(int), 1, temp_mapping_output_file);
@@ -1253,7 +1256,7 @@ inline void SAMOutputTools<SAMMapping>::AppendMapping(uint32_t rid, const Sequen 
  //std::string strand = (mapping.direction & 1) == 1 ? "+" : "-";

  //uint32_t mapping_end_position = mapping.fragment_start_position + mapping.fragment_length;

  const char *reference_sequence_name = (mapping.flag & BAM_FUNMAP) > 0 ? "*" : reference.GetSequenceNameAt(rid);

  this->AppendMappingOutput(mapping.read_name + "\t" + std::to_string(mapping.flag) + "\t" + std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping.mapq) + "\t" + mapping.GenerateCigarString() + "\t*\t" + std::to_string(0) + "\t" + std::to_string(0) + "\t*\t*\t" + mapping.GenerateIntTagString("NM", mapping.NM) + "\tMD:Z:" + mapping.MD + "\n");

  this->AppendMappingOutput(mapping.read_name + "\t" + std::to_string(mapping.flag) + "\t" + std::string(reference_sequence_name) + "\t" + std::to_string(mapping.GetStartPosition()) + "\t" + std::to_string(mapping.mapq) + "\t" + mapping.GenerateCigarString() + "\t*\t" + std::to_string(0) + "\t" + std::to_string(0) + "\t*\t*\t" + mapping.GenerateIntTagString("NM", mapping.NM) + "\tMD:Z:" + mapping.MD + "\tCB:Z:" + Seed2Sequence(mapping.cell_barcode, cell_barcode_length_) + "\n");

}



template <>

CRA Git | Maintained and supported by SUSTech CRA and CCSE