Grcov report - read

1

#[doc(inline)]

2

use std::path::PathBuf;

3

4

use color_eyre::{eyre::eyre, Result};

5

use ndarray::{Array2, ShapeBuilder};

6

use rust_htslib::bcf::header::HeaderView;

7

use rust_htslib::bcf::record::GenotypeAllele;

8

use rust_htslib::bcf::IndexedReader;

9

use rust_htslib::bcf::Read;

10

use rust_htslib::bcf::Record;

11

12

use crate::args::StandardArgs;

13

use crate::io::read_sample_ids;

14

use crate::structs::{Coord, PhasedMatrix};

15

use crate::utils::filter_samples;

16

17

888

pub fn get_reader(

18

888

    path: &PathBuf,

19

888

    contig: &str,

20

888

    range: Option<(u64, u64)>,

21

888

) -> Result<IndexedReader> {

22

888

    let mut reader = IndexedReader::from_path(path)?;

23

888

    let rid = reader.header().name2rid(contig.as_bytes())?;

24

25

888

    match range {

26

        // RUST-HTSLIB is 0-based so subtract 1

27

54

        Some((start, end)) => reader.fetch(rid, start - 1, Some(end - 1))?,

28

834

        None => reader.fetch(rid, 0, None)?,

29

};

30

31

888

    Ok(reader)

32

888

33

34

870

pub fn get_samples(header: &HeaderView) -> Result<Vec<String>> {

35

870

    header

36

870

        .samples()

37

870

        .into_iter()

38

12168

        .map(|sample| Ok(String::from_utf8_lossy(sample).to_string()))

39

870

        .collect()

40

870

41

42

/// Creates a two dimensional array that contains samples and variants. REF homozygotes are denoted

43

/// by 0.  Heterozygotes are denoted by 1. ALT homozygotes are denoted by 2.

44

18

pub fn read_vcf_to_contra_matrix(

45

18

    args: &StandardArgs,

46

18

    contig: &str,

47

18

    range: Option<(u64, u64)>,

48

18

) -> Result<(Array2<u8>, Vec<i64>)> {

49

18

    let mut reader = get_reader(&args.file, contig, range)?;

50

18

    let samples = get_samples(reader.header())?;

51

18

    let wanted = read_sample_ids(&args.samples)?;

52

18

    let sample_indexes = filter_samples(&samples, wanted);

53

18

54

18

    let mut pos_coords = vec![];

55

18

    let mut haplotypes = vec![];

56

18

    let mut records_n = 0;

57

58

804

    for (_i, record_result) in reader.records().enumerate() {

59

804

        let record = record_result?;

60

61

804

        if record.alleles().len() > 2 {

62

6

            return Err(eyre!("VCF has multiallelics. Use bcftools norm."));

63

798

64

798

65

798

        // RUST-HTSLIB is 0-based so add 1

66

798

        let pos = record.pos() + 1;

67

798

        pos_coords.push(pos);

68

798

69

798

        if check_info_limit(args, &record, contig, pos)? {

70

            // Loop genotypes and assign to matrix as 0, 1 or 2

71

798

            let gts = record.genotypes()?;

72

73

11970

            for i in &sample_indexes {

74

22344

                let gt_vec: Vec<Option<u32>> = gts.get(*i).iter().map(|gt| gt.index()).collect();

75

11172

                if gt_vec.len() > 2 {

76

                    panic!("Three genotypes present at {pos}");

77

11172

78

11172

                match (gt_vec[0], gt_vec[1]) {

79

10302

                    (Some(0), Some(0)) => haplotypes.push(0),

80

                    (Some(1), Some(1)) => haplotypes.push(2),

81

                    (Some(2), Some(_)) | (Some(_), Some(2)) => normalize_error(pos)?,

82

870

                    _ => haplotypes.push(1),

83

84

85

798

            records_n += 1;

86

87

88

89

12

    let matrix = Array2::from_shape_vec((sample_indexes.len(), records_n).f(), haplotypes)?;

90

91

12

    Ok((matrix, pos_coords))

92

18

93

94

/// Read a tabix indexed VCF using the contig name given. Create a 2-D genotype matrix where both

95

/// alleles are stacked in the same matrix. Rows are sample alleles and columns are variants. To

96

/// find the second genotype of the same sample use n + vec.len()

97

/// If a variant presents any missing genotypes, the entire variant is ignored

98

///

99

228

pub fn read_vcf_to_matrix(

100

228

    args: &StandardArgs,

101

228

    contig: &str,

102

228

    variant_pos: i64,

103

228

    range: Option<(u64, u64)>,

104

228

) -> Result<PhasedMatrix> {

105

228

    let mut reader = get_reader(&args.file, contig, range)?;

106

228

    let samples = get_samples(reader.header())?;

107

108

228

    let wanted = read_sample_ids(&args.samples)?;

109

228

    let sample_indexes = filter_samples(&samples, wanted);

110

228

111

228

    let mut markers_0 = vec![];

112

228

    let mut markers_1 = vec![];

113

228

114

228

    let mut coords = vec![];

115

228

116

228

    let mut records_n = 0;

117

228

    let mut variant_idx = 0;

118

228

119

228

    let mut gt_buffer = rust_htslib::bcf::record::Buffer::new();

120

228

    tracing::debug!("Reading contig {contig} at position {variant_pos}");

121

122

12240

    for record in reader.records() {

123

12240

        let record = record?;

124

        // RUST-HTSLIB is 0-based so add 1

125

12240

        let pos = record.pos() + 1;

126

12240

127

12240

        tracing::trace!("Reading record at position {pos}");

128

12240

        if pos >= variant_pos && variant_idx == 0 {

129

234

            variant_idx = records_n;

130

12018

131

132

12252

        if check_info_limit(args, &record, contig, pos)? {

133

11862

            let gts = record.genotypes_shared_buffer(&mut gt_buffer)?;

134

11862

            let mut phased_counter = 0;

135

11862

            let mut missing = false;

136

11862

            let mut pending_for_0 = vec![];

137

11862

            let mut pending_for_1 = vec![];

138

139

177978

            for i in &sample_indexes {

140

332700

                for (ii, gta) in gts.get(*i).iter().enumerate() {

141

332700

                    match gta {

142

165906

                        GenotypeAllele::Unphased(v) => {

143

165906

                            if v == &2 {

144

6

                                normalize_error(pos)?;

145

165900

146

147

165900

                            match ii {

148

165900

                                0 => pending_for_0.push(*v as u8),

149

                                1 => pending_for_1.push(*v as u8),

150

                                _ => panic!("three genotypes in a variant at {pos}"),

151

152

153

166794

                        GenotypeAllele::Phased(v) => {

154

166794

                            if v == &2 {

155

                                normalize_error(pos)?;

156

166794

157

166794

                            phased_counter += 1;

158

166794

                            match ii {

159

408

                                0 => pending_for_0.push(*v as u8),

160

166386

                                1 => pending_for_1.push(*v as u8),

161

                                _ => panic!("three genotypes in a variant at {pos}"),

162

163

164

                        GenotypeAllele::UnphasedMissing => missing = true,

165

                        GenotypeAllele::PhasedMissing => missing = true,

166

167

168

169

11922

            if !missing {

170

11922

                records_n += 1;

171

11922

                markers_0.extend(pending_for_0);

172

11922

                markers_1.extend(pending_for_1);

173

11922

174

11922

                let alleles = record.alleles();

175

11922

                let reference = String::from_utf8_lossy(alleles.get(0).unwrap()).to_string();

176

11922

                let alt = String::from_utf8_lossy(alleles.get(1).unwrap()).to_string();

177

11922

178

11922

                let coord = Coord {

179

11922

                    contig: contig.to_string(),

180

11922

                    reference,

181

11922

                    alt,

182

11922

                    pos,

183

11922

                    phase: phased_counter == samples.len(),

184

11922

};

185

11922

                coords.push(coord);

186

11922

            } else {

187

                tracing::warn!("Genotypes missing at {contig}:{pos}, these variants are not accounted in any analyses");

188

189

384

190

191

192

222

    if variant_idx == 0 {

193

        return Err(eyre!(

194

            "The given coordinate is larger than the largest found position. Total records read: {records_n}. Check your coordinates and vcf file. Comparing to a haplotype also limits min and max coordinates."

195

));

196

222

197

198

222

    let matrix_0 = Array2::from_shape_vec((sample_indexes.len(), records_n).f(), markers_0)?;

199

222

    let matrix_1 = Array2::from_shape_vec((sample_indexes.len(), records_n).f(), markers_1)?;

200

201

222

    let matrix = ndarray::concatenate!(ndarray::Axis(0), matrix_0, matrix_1);

202

222

203

222

    let samples: Vec<String> = sample_indexes

204

222

        .iter()

205

222

        .cycle()

206

6216

        .map(|s| &samples[*s])

207

222

        .take(sample_indexes.len() * 2)

208

222

        .cloned()

209

222

        .collect();

210

222

211

222

    Ok(PhasedMatrix::new(variant_idx, matrix, samples, coords))

212

228

213

214

/// Read a tabix indexed VCF using the contig name given. Create a 2-D genotype matrix where both

215

/// alleles are stacked in the same matrix. Rows are sample alleles and columns are variants. To

216

/// find the second genotype of the same sample use n + vec.len()

217

/// Denote missing genotypes as 99

218

12

pub fn read_vcf_to_matrix_missing_allowed(

219

12

    args: &StandardArgs,

220

12

    contig: &str,

221

12

    variant_pos: i64,

222

12

    range: Option<(u64, u64)>,

223

12

) -> Result<PhasedMatrix> {

224

12

    let mut reader = get_reader(&args.file, contig, range)?;

225

12

    let samples = get_samples(reader.header())?;

226

227

12

    let wanted = read_sample_ids(&args.samples)?;

228

12

    let sample_indexes = filter_samples(&samples, wanted);

229

12

230

12

    let mut markers_0 = vec![];

231

12

    let mut markers_1 = vec![];

232

12

233

12

    let mut coords = vec![];

234

12

235

12

    let mut records_n = 0;

236

12

    let mut variant_idx = 0;

237

12

238

12

    let mut gt_buffer = rust_htslib::bcf::record::Buffer::new();

239

240

426

    for record in reader.records() {

241

426

        let record = record?;

242

        // RUST-HTSLIB is 0-based so add 1

243

426

        let pos = record.pos() + 1;

244

426

245

426

        if pos >= variant_pos && variant_idx == 0 {

246

6

            variant_idx = records_n;

247

420

248

249

426

        if check_info_limit(args, &record, contig, pos)? {

250

426

            let gts = record.genotypes_shared_buffer(&mut gt_buffer)?;

251

426

            let mut phased_counter = 0;

252

253

6324

            for i in &sample_indexes {

254

11802

                for (ii, gta) in gts.get(*i).iter().enumerate() {

255

11802

                    match gta {

256

5892

                        GenotypeAllele::Unphased(v) => {

257

5892

                            if v == &2 {

258

6

                                normalize_error(pos)?;

259

5886

260

5886

                            match ii {

261

5886

                                0 => markers_0.push(*v as u8),

262

                                1 => markers_1.push(*v as u8),

263

                                _ => panic!("three genotypes in a variant at {pos}"),

264

265

266

5892

                        GenotypeAllele::Phased(v) => {

267

5892

                            if v == &2 {

268

                                normalize_error(pos)?;

269

5892

270

5892

                            phased_counter += 1;

271

5892

                            match ii {

272

                                0 => markers_0.push(*v as u8),

273

5892

                                1 => markers_1.push(*v as u8),

274

                                _ => panic!("three genotypes in a variant at {pos}"),

275

276

277

12

                        GenotypeAllele::UnphasedMissing => match ii {

278

12

                            0 => markers_0.push(99),

279

                            1 => markers_1.push(99),

280

                            _ => panic!("three genotypes in a variant at {pos}"),

281

},

282

6

                        GenotypeAllele::PhasedMissing => match ii {

283

                            0 => markers_0.push(99),

284

6

                            1 => markers_1.push(99),

285

                            _ => panic!("three genotypes in a variant at {pos}"),

286

},

287

288

289

290

420

            records_n += 1;

291

420

292

420

            let alleles = record.alleles();

293

420

            let reference = String::from_utf8_lossy(alleles.get(0).unwrap()).to_string();

294

420

            let alt = String::from_utf8_lossy(alleles.get(1).unwrap()).to_string();

295

420

296

420

            let coord = Coord {

297

420

                contig: contig.to_string(),

298

420

                reference,

299

420

                alt,

300

420

                pos,

301

420

                phase: phased_counter == samples.len(),

302

420

};

303

420

            coords.push(coord);

304

305

306

307

6

    if variant_idx == 0 {

308

        return Err(eyre!(

309

            "The given coordinate is larger than the largest found position. Total records read: {records_n}. Check your coordinates and vcf file. Comparing to a haplotype also limits min and max coordinates."

310

));

311

6

312

313

6

    let matrix_0 = Array2::from_shape_vec((sample_indexes.len(), records_n).f(), markers_0)?;

314

6

    let matrix_1 = Array2::from_shape_vec((sample_indexes.len(), records_n).f(), markers_1)?;

315

316

6

    let matrix = ndarray::concatenate!(ndarray::Axis(0), matrix_0, matrix_1);

317

6

318

6

    let samples: Vec<String> = sample_indexes

319

6

        .iter()

320

6

        .cycle()

321

168

        .map(|s| &samples[*s])

322

6

        .take(sample_indexes.len() * 2)

323

6

        .cloned()

324

6

        .collect();

325

6

326

6

    Ok(PhasedMatrix::new(variant_idx, matrix, samples, coords))

327

12

328

329

pub fn check_info_limit(

330

    args: &StandardArgs,

331

    record: &Record,

332

    contig: &str,

333

    pos: i64,

334

) -> Result<bool> {

335

13464

    if let Some(info_limit) = args.info_limit {

336

744

        match record.info(b"INFO").float()? {

337

756

            Some(buffer) => Ok(buffer[0] >= info_limit),

338

            None => return Err(eyre!("No INFO score on a variant at: {contig}:{pos}")),

339

340

    } else {

341

12720

        Ok(true)

342

343

13476

344

345

12

fn normalize_error(pos: i64) -> Result<()> {

346

12

    return Err(eyre!(

347

12

        "At pos {pos} allele count != 2. Normalize alleles using bcftools norm"

348

12

));

349

12