-
Notifications
You must be signed in to change notification settings - Fork 1
/
create_vcf_file.R
32 lines (25 loc) · 1.13 KB
/
create_vcf_file.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
library(data.table)
library(tidyverse)
haps <- fread(snakemake@input$haplotypes)
#genomes_1000 <- fread(snakemake@config$thousand_genomes_snps)
#names(genomes_1000) <- c("chromosome", "position", "ref", "alt")
#haps <- merge(haps, genomes_1000, on = c("chromosome", "position"), all.x = T)
haps <- select(haps, chromosome, position, ref, alt) %>%
filter(alt != "-")
#rename columns and add columns required by VCF
haps$`#CHROM` <- haps$chromosome
haps$POS <- haps$position
haps$REF <- haps$ref
haps$ALT <- haps$alt
haps$QUAL <- "."
haps$ID <- "."
haps$FILTER <- "PASS"
haps$INFO <- "."
haps <- select(haps, `#CHROM`, POS, ID, REF, ALT,QUAL,FILTER,INFO) %>%
distinct()
write("##fileformat=VCFv4.2", snakemake@output$haplotypesvcf)
fwrite(haps, file = snakemake@output$haplotypesvcf, sep = "\t", append = T, col.names = T)
# fread("../../../schnapps_input_pipeline/metadata/samples.csv") %>%
# .[, cell_id := str_remove(cell_id, "/work/shah/users/william1/projects//schnapps_input_pipeline/testdata//")] %>%
# .[, bamfiles := str_replace(bamfiles, "\\//", "/")] %>%
# fwrite("../../../schnapps_input_pipeline/metadata/samples.csv")