Make chromosome plots

[adamjorr]

Make chromosome plots and put mutations on plot.

[adamjorr]

To get gff of genes on 1st 11 chr:

zcat ~/Downloads/Egrandis_201_v1.1.gene.gff3.gz | sed -n -E '/^scaffold_[0-9]\t|scaffold_1[0,1]\t/p' | awk '$3 == "gene" {print}' | less -S

[adamjorr]

zcat ~/Downloads/Egrandis_201_v1.1.gene.gff3.gz | sed -n -E '/^scaffold_[0-9]\t|scaffold_1[0,1]\t/p' | awk '$3 == "gene" {print}' | cut -f1,4,5,9 | awk 'BEGIN{OFS="\t";print "chr","start","end","name","gieStain"}{gsub(/.+[nN]ame=/,"",$4); print $0, "gpos50"}' | less -S

[adamjorr]

To get lengths of first 11 chr:

head -n 11 ref.fa.fai | cut -f1,2 | awk 'BEGIN{OFS="\t";print "chr","start","end"}{print $1
,"1",$2}' | less -S

https://bernatgel.github.io/karyoploter_tutorial//Tutorial/CustomGenomes/CustomGenomes.html

[adamjorr]

to get lengths from first 11 chr from the header of a vcf:

bcftools view -h filtered_bed_excluded.vcf | sed -n -E '/ID=scaffold_[0-9],|ID=scaffold_1[01],/p' | tr -d '<>' | sed -nE 's/length=//gp' | cut -d= -f3 | tr , '\t' | awk 'BEGIN{OFS="\t";print "chr","start","end"}{print $1, 1, $2}' | less -S

[adamjorr]

library(karyoploteR)
library(tidyverse)

genome <- toGRanges("chromosomes.txt")
bands <- toGRanges("genes.txt")
withrepeats <- read_tsv("positions_with_repeats.txt", col_names=c("chr","pos"))
norepeats <- read_tsv("positions_no_repeats.txt", col_names=c("chr","pos"))

pdf("chromosome_plot.pdf",10,8)
kp <- plotKaryotype(genome = genome, cytobands = bands, cex=.75)

# plot_above <- function(chr, pos, kp, y, ...){
#   kpPoints(kp, chr=chr, x=pos, y=y, ...)
# }
# norepeats %>%
#   pwalk(plot_above, kp=kp, y = .5, col = "green")
# 
# withrepeats %>%
#   pwalk(plot_above, kp=kp, y = .2)

kpPoints(kp, chr=withrepeats$chr, x=withrepeats$pos, y=.2)
kpPoints(kp, chr=norepeats$chr, x=norepeats$pos, y = .5, col = "green")
dev.off()

variant_table <- norepeats %>%
  mutate(repeatregion = FALSE) %>%
  right_join(withrepeats) %>%
  mutate(repeatregion = is.na(.$repeatregion))

write_tsv(variant_table, "variant_table.tsv")

[adamjorr]

chromosome_plot.pdf

[adamjorr]

library(karyoploteR)
library(VariantAnnotation)
library(IRanges)
library(GenomicRanges)
library(tidyverse)


genome <- toGRanges("chromosomes.txt")
bands <- toGRanges("genes.txt")
withrepeats <- read_tsv("positions_with_repeats.txt", col_names=c("chr","pos"))
norepeats <- read_tsv("positions_no_repeats.txt", col_names=c("chr","pos"))
genes <- read_tsv("genes.txt") %>%
  group_by(chr)

pdf("chromosome_plot.pdf",10,8)
kp <- plotKaryotype(genome = genome, cytobands = bands, cex=.75)

# plot_above <- function(chr, pos, kp, y, ...){
#   kpPoints(kp, chr=chr, x=pos, y=y, ...)
# }
# norepeats %>%
#   pwalk(plot_above, kp=kp, y = .5, col = "green")
# 
# withrepeats %>%
#   pwalk(plot_above, kp=kp, y = .2)

kpPoints(kp, chr=withrepeats$chr, x=withrepeats$pos, y=.2)
kpPoints(kp, chr=norepeats$chr, x=norepeats$pos, y = .5, col = "green")
dev.off()

in_range_ <- function(chrm, pos, genes){
  rnges <- genes %>% filter(chr == chrm)
  btwn <- function(pos, start, end){between(pos,start,end)}
  result <- map2(rnges$start,rnges$end, ~ btwn(pos,.x,.y))
  return(TRUE %in% result)
}
in_range <- function(chrm,pos,genes){
  return(map2(chrm,pos,~ in_range_(.x,.y,genes)))
}

variant_table <- norepeats %>%
  mutate(repeatregion = FALSE) %>%
  right_join(withrepeats) %>%
  mutate(repeatregion = is.na(.$repeatregion)) %>%
  mutate(ingene = unlist(in_range(chr,pos,genes)))

write_tsv(variant_table, "variant_table.tsv")

# this may break if tidyverse is loaded??
vcf <- readVcf("passed.vcf","e_mel_3.fa")
gff <- GenomicFeatures::makeTxDbFromGFF("e_mel_3_genes.gff3")
locs <- locateVariants(vcf,gff,AllVariants())
fa <- FaFile("e_mel_3.fa")
indexFa("e_mel_3.fa")
coding <- predictCoding(vcf,gff,seqSource=fa)

# load tidy things here
positions <- as_tibble(start(ranges(coding)))
names(positions) <- "pos"
chrs <- as_tibble(seqnames(coding))
names(chrs) <- "chr"
coding <- chrs %>%
  bind_cols(positions) %>%
  bind_cols(as_tibble(mcols(coding)))

positions <- as_tibble(start(ranges(locs)))
names(positions) <- "pos"
chrs <- as_tibble(seqnames(locs))
names(chrs) <- "chr"

locs <- chrs %>%
  bind_cols(positions) %>%
  bind_cols(as_tibble(mcols(locs)))

type <- locs %>%
  select(chr, pos, LOCATION) %>%
  distinct()

coding <- coding %>%
  # mutate(ALT = unlist(map(coding$ALT, ~ as.character(.))))
  select(chr, pos, CONSEQUENCE, REFAA, VARAA)

variant_table <- variant_table %>%
  left_join(type,by = c("chr","pos")) %>%
  left_join(coding,by = c("chr","pos")) 

# variant_table <- variant_table %>%
#   mutate(seq = unlist(map2(.$chr,.$pos, ~ as.character(getSeq(fa, GRanges(seqnames = .x, IRanges(start = .y - 1000, end = .y + 1000)))))))

gts <- as.tibble(geno(genotypeCodesToNucleotides(vcf))$GT)

refalt <- as_tibble(mcols(rowRanges(vcf))) %>%
  select(REF, ALT) %>%
  mutate(ALT = unlist(map(.$ALT, ~ as.character(.)))) %>%
  bind_cols(as_tibble(as.data.frame(seqnames(rowRanges(vcf))))) %>%
  rename(chr = value) %>%
  bind_cols(as_tibble(start(ranges(rowRanges(vcf))))) %>%
  rename(pos = value) %>%
  select(chr, pos, REF, ALT) %>%
  bind_cols(gts) %>%
  mutate(upstreamkb = unlist(map2(.$chr,.$pos, ~as.character(getSeq(fa,GRanges(seqnames = .x, IRanges(start = .y - 1000, end = .y - 1))))))) %>%
  mutate(downstreamkb = unlist(map2(.$chr,.$pos, ~as.character(getSeq(fa,GRanges(seqnames = .x, IRanges(start = .y + 1 , end = .y + 1000))))))) %>%
  select(chr:ALT, upstreamkb, M1a:M8c, downstreamkb)

variant_table <- variant_table %>%
  left_join(refalt, by=c("chr","pos"))

write_tsv(variant_table,"variants.tsv")