commands.txt

kinit 
klist 

pwd
ls

cd /filer/projekte/gbs_course 
cd data
ls
ls -lhS gbs_reads | less
man ls
cd gbs_reads
ls | wc -l
zcat Morex.fastq.gz | less -S
cd ..
du -ch gbs_reads

less MxB_92_RILs_phenotypes.csv

cd /filer/projekte/gbs_course
mkdir -p mascher/gbs
cd mascher/gbs

pwd

mkdir gbs_reads
datadir="/filer/projekte/gbs_course/data/"
ln -s $datadir/gbs_reads/*fastq.gz gbs_reads
ls -l gbs_reads | less -S
ls gbs_reads | wc -l

ln -s $datadir/MxB_92_RILs_phenotypes.csv

zcat gbs_reads/Morex.fastq.gz | wc -l > Morex.len

for i in gbs_reads/*fastq.gz; do
 name=`basename $i | cut -d . -f 1`
 count=`zcat $i | wc -l`
 echo $name $count
done > fastq_line_counts.txt

cat fastq_line_counts.txt | awk '{print $1,$2/4}' | sort -nk 2 > fastq_read_counts.txt

module load bwa samtools
module list

ref='/filer/projekte/gbs_course/data/reference/Hordeum_vulgare_MorexV3_assembly.fasta'
bwa shm $ref
bwa shm -l

zcat gbs_reads/Morex.fastq.gz | head -n 40000 > Morex_10k.fastq
bwa mem $ref Morex_10k.fastq > Morex_10k.sam 
samtools sort Morex_10k.sam > Morex_10k.bam 
samtools view Morex_10k.bam | less -S

bwa mem $ref Morex_10k.fastq | samtools sort > Morex_10k.bam

k5reauth tmux

tmux ls
tmux rename mapping
tmux ls
tmux detach 
tmux ls
tmux attach -t mapping

ref='/filer/projekte/gbs_course/data/reference/Hordeum_vulgare_MorexV3_assembly.fasta'
for i in gbs_reads/*fastq.gz; do
 name=`echo $i | cut -d . -f 1`
 bwa mem -t 4 $ref $i | samtools sort > $name.bam
done 2> bwa.err

find gbs_reads | grep -c 'bam$'
top -u mascher

for i in gbs_reads/*bam; do
 name=`basename $i | cut -d . -f 1`
 count=`samtools view -cq 20 $i`
 echo $name $count
done > mapped_reads.txt

LC_ALL=C sort fastq_read_counts.txt > tmp1 
LC_ALL=C sort mapped_reads.txt > tmp2 

join tmp1 tmp2 | awk '{print $0,$3/$2*100}' | sort -nk 4 | column -t > mapping_stats.txt // <2>

rm -f tmp1 tmp2

column -t mapping_stats.txt  | less -S

ls gbs_reads/*bam | sort > bam_list.txt

tmux
module load bcftools
ref='/filer/projekte/gbs_course/data/reference/Hordeum_vulgare_MorexV3_assembly.fasta'
#bamlist='bam_list.txt' 
bamlist='/filer-dg/agruppen/dg6/mascher/DG/mxb_course_221201/try_221216/bam_list_Martin.txt' 
vcf='bcftools_SNP_calling.vcf'

bcftools mpileup --skip-indels --fasta-ref $ref --bam-list $bamlist --min-MQ 20 --annotate AD \
 | bcftools view -i 'INFO/DP > 100' \
 | bcftools call --consensus-caller --variants-only --output $vcf

filter='/filer/projekte/gbs_course/scripts/filter_vcf.zsh'
#vcf='bcftools_SNP_calling.vcf'
vcf='/filer-dg/agruppen/dg6/mascher/DG/mxb_course_221201/try_221216/bcftools_SNP_calling.vcf' 
fvcf='bcftools_SNP_calling_filtered.vcf'

$filter --vcf $vcf --dphom 2 --dphet 4 --minmaf 0.2 --minpresent 0.9 --minhomp 0.9 > $fvcf

less MxB_92_RILs_phenotypes.csv
bcftools query -l $fvcf | cut -d / -f 2 | cut -d . -f 1 > new_sample_names.txt
bcftools reheader -s new_sample_names.txt $fvcf > bcftools_SNP_calling_filtered_newNames.vcf
grep -v '^##' bcftools_SNP_calling_filtered_newNames.vcf | column -t | less -S

module load R/4.1.1
R 

.libPaths("/filer/projekte/gbs_course/Rlibs/4.1.1")

library("qtl")
library("ASMap")
library("utl")

vcf <- 'bcftools_SNP_calling_filtered_newNames.vcf'
ids <- read.table("new_sample_names.txt", head=F)[, 1]
genfile <- 'bcftools_SNP_calling_geno.csv'
founders <- c("Morex", "Barke")
samples <- setdiff(ids, founders)

convert_vcf_to_genfile(vcf, genfile, samples, founders)

cmd <- "sed -Ei '1{p; s/id|chr.H://g; h; d}; 2G'"
paste(cmd, genfile) |> system()

pheno <- 'MxB_92_RILs_phenotypes.csv'
read.cross(format="csvs", genfile=genfile, phefile=pheno, crosstype= "f2",  genotypes=c("AA","AB","BB")) -> mxb
convert2bcsft(mxb, BC.gen=0, F.gen=8,  estimate.map=F) -> mxb
summary(mxb)
summary.map(mxb)

nind(mxb) # number of inviduals
nchr(mxb) # number of chromosomes (actually linkage group)
totmar(mxb) # total number of markers
nmar(mxb) # number of markers per chromosome
nphe(mxb) # number of phenotypes

pdf("plot_pheno.pdf")
plotPheno(mxb, pheno.col=1, xlab="phenotypes")
plotPheno(mxb, pheno.col=2)
dev.off()

rescalemap(mxb, 1/1e6) -> mxb
summary.map(mxb)

pdf("plot_map.pdf")
plot.map(mxb, main="Physical map", ylab="Location (Mb)")
dev.off()

copy(mxb) -> mxb_physical
saveRDS(mxb_physical, "mxb_physical.Rds") 
readRDS("mxb_physical.Rds") -> mxb_physical

summary.map(mstmap(mxb, id="id"))

findDupMarkers(mxb, exact.only=FALSE) -> dups
unlist(dups) |> length()

mxb <- drop.markers(mxb, unlist(dups))
summary(mxb)

cg <- comparegeno(mxb)

pdf("compare_geno.pdf")
hist(cg[lower.tri(cg)], breaks=seq(0, 1, len=101), xlab="No. matching genotypes")
dev.off()

summary(cg, thresh=0.9) -> dups
subset(mxb, ind=setdiff(samples, dups$inds2)) -> mxb

pdf("plot.geno.pdf")
geno.image(mxb, col=c("white", "red", "violet", "blue"))
dev.off()

ntyped(mxb) |> sort() |> head(n=1) -> rm.ind

subset(mxb, ind=setdiff(names(ntyped(mxb)), names(rm.ind))) -> mxb
summary(mxb)

set.seed(1)
lapply(pull.map(mxb), function(i) pickMarkerSubset(i, 1)) |> unlist() -> keep

drop.markers(mxb, setdiff(markernames(mxb), keep))  -> mxb

pdf("plot.geno_1Mb.pdf")
geno.image(mxb, col=c("white", "red", "violet", "blue"))
dev.off()

mstmap.cross(mxb, id="id") -> mxb
summary.map(mxb)

pdf("align_maps.pdf")
alignCross(mxb, maps=list(mxb_physical), layout=c(1,1))
dev.off()

pull.map(mxb, as.table=T) -> a
pull.map(mxb_physical, as.table=T) -> b
merge(a, b, by=0) -> m 
sapply(split(m, m$chr.x), function(i) with(i, cor(pos.x, pos.y, method='s'))) -> cc

names(which(cc < 0)) -> flip.chr
flip.order(mxb, flip.chr) -> mxb

pdf("align_maps_2.pdf")
alignCross(mxb, maps=list(mxb_physical), layout=c(1,1))
dev.off()

pdf("plot.geno.mstmap.pdf")
geno.image(mxb, col=c("white", "red", "violet", "blue"))
dev.off()

rescalemap(mxb_physical, 1/5) -> mxb_physical

pdf("plot_map_2.pdf")
plotMap(mxb, mxb_physical)
dev.off()

copy(mxb) -> mxb_genetic
saveRDS(mxb_genetic, file="mxb_genetic.Rds")

.libPaths("/filer/projekte/gbs_course/Rlibs/4.1.1") // <1>

library("qtl") 
readRDS("mxb_genetic.Rds") -> mxb

calc.genoprob(mxb) -> mxb

mxb$pheno$row_type <- ifelse(mxb$pheno$row_type == 2, 0, 1)

scanone(mxb, pheno.col=1, method="mr", model="binary") -> out
summary(out)

pdf("plot_row_type.pdf")
plot(out)
dev.off()

operm <- scanone(mxb, method="mr", n.perm=1000, pheno.col=1) 

summary(operm, alpha=c(0.05, 0.2))

pdf("plot_operm.pdf")
hist(as.numeric(operm))
abline(v=summary(operm, alpha=c(0.05)), col='red')
dev.off()

summary(out, perms=operm, pvalues=TRUE)

pdf("plot_row_type_threshold.pdf")
plot(out)
abline(h=3.5, col='red')
dev.off()

lodint(out, chr="2H", drop=2)

rownames(max(out)) -> mar
pdf("plot_pxg.pdf")
plotPXG(mxb, pheno.col=1, marker=mar)
dev.off()

.libPaths("/filer/projekte/gbs_course/Rlibs/4.1.1") 
library(data.table)
library(xlsx)
library("qtl") 

fread('MxB_92_RILs_phenotypes.csv')->p

data.table(read.xlsx('/filer-dg/agruppen/dg6/mascher/DG/mxb_course_221201/prep/KASP_results_PGR23.xlsx', 1))->k
setnames(k, "Samples", "id")
k[grepl("F11$", id)][, id := sub("MBRIL-", "07-", sub("_F11$", "", id))][] -> k

k[p, on="id"] -> m
m[, .(row_type, plant_height, RAW1_Henry, VRS1_Henry, RAW1_ZsaZsa, VRS1_ZsaZsa, RAW1_Tania, VRS1_Tania, id)] -> m

unique(m[, .(row_type, VRS1_Tania)])

readRDS('mxb_genetic.Rds')

data.frame(m[id %in% mxb$pheno$id]) -> mxb$pheno
mxb
colnames(mxb$pheno)

pdf("plot_KASP.pdf")
plotPheno(mxb, pheno.col=3)
plotPheno(mxb, pheno.col=4)
plotPheno(mxb, pheno.col=5)
plotPheno(mxb, pheno.col=6)
plotPheno(mxb, pheno.col=7)
plotPheno(mxb, pheno.col=8)
dev.off()

mxb$pheno$VRS1_Tania <- ifelse(mxb$pheno$VRS1_Tania == "C", 0, 1)
scanone(mxb, pheno.col=8, method="mr", model="binary") -> out

pdf("plot_VRS1_Tania.pdf")
plot(out)
dev.off()

quit()