#########################################################
# AprendeBioinformaticaEncasa Parte V
# Herramientass de linea de comnnados UNIX
# para la manipulación y análisis de datos genomicos (II)
#########################################################

#Navegar a nuestro directorio
cd /mnt/c/.../mequedoencasa

# Visualizar archivo BAM
samtools view bam/1M68_pH5_0.04CO2_R1.bam | less

#Ver estadísticas de BAM
samtools flagstats bam/1M68_pH5_0.04CO2_R1.bam

# Contar lecturas alineadas con un estado de emparejamiento "correcto"
samtools view -f 2 bam/1M68_pH5_0.04CO2_R1.bam | wc -l

# Eliminar duplicados 
samtools sort -n -o bam/1M68_pH5_0.04CO2_R1.namesort.bam bam/1M68_pH5_0.04CO2_R1.bam
samtools fixmate -m bam/1M68_pH5_0.04CO2_R1.namesort.bam bam/1M68_pH5_0.04CO2_R1.namesort.fixmate.bam
samtools sort -o bam/1M68_pH5_0.04CO2_R1.namesort.fixmate.sorted.bam bam/1M68_pH5_0.04CO2_R1.namesort.fixmate.bam

# Eliminar duplicados concatenando comandos para evitar archivos intermedios
samtools sort -n bam/1M68_pH5_0.04CO2_R1.bam |  samtools fixmate -m - - | samtools sort - | samtools markdup - bam/1M68_ph5_0.04CO2_R1.namesort.fixmate.sorted.markdup.bam

# Interseccion de regiones en dos archivos BED
bedtools intersect -a bed/regions_example.bed -b bed/regions_example2.bed > bed/intersect_1_2.bed

# Interseccion de regiones de un archivo BED con otros dos archivos BED
bedtools intersect -a bed/regions_example.bed -b bed/regions_example2.bed bed/regions_example3.bed > bed/intersect_1_2-3.bed

# Contar número de intersecciones de un archivo BED con regiones de un segundo archivo
bedtools intersect -c -a bed/regions_example.bed -b bed/regions_example2.bed > bed/regions_exampleIn2.bed

# Unir regiones de varios archivos BED (Debemos ordenar las regiones)
cat bed/regions_example.bed bed/regions_example2.bed  bed/regions_example3.bed  > bed/all_regions.bed
bedtools sort bed/all_regions.bed > bed/all_regions.sorted.bed
sort -k1,1 -k2,2n bed/all_regions.bed > bed/all_regions.sorted.bed

# Juntar regiones solapantes en un archivo BED 
bedtools merge -i bed/all_regions.sorted.bed > bed/all_regions_merged.bed

# Juntar regiones solapantes (o que estén a menos de 20nt) en un archivo BED 
bedtools merge -d 20 -i bed/all_regions.sorted.bed > bed/all_regions_mergedAt20nt.bed

# Juntar regiones solapantes en un archivo BED teniendo en cuenta la direccionalidad
bedtools merge -s -i bed/all_regions.sorted.bed > bed/all_regions_merged_samestrand.bed

# Crear archivo de tamaño de cromosoma
grep -P '\tchromosome\t' genes/Saccharomyces_cerevisiae.R64-1-1.99.gff3 | cut -f1,5 > genes/Saccharomyces_cerevisiae.chromSizes

# Crear archivo de regiones 500bp antes de cada gen
bedtools flank -i bed/s_cerevisiae_genes.bed -g genes/Saccharomyces_cerevisiae.chromSizes -s -r 0 -l 500 > bed/s_cerevisiae_upstream500.bed

# Calcular coverage para cada elemento de un archivo BED de un único alineamiento
mkdir bedcoverage
bedtools coverage -a bed/s_cerevisiae_genes.bed -b bam/1M68_pH5_0.04CO2_R1.bam > bedcoverage/genes_1M68_pH5_0.04CO2_R1.tsv

# Obtener coverage para cada elemento de un archivo BED de múltiples archivos de alineamiento
mkdir bedcoverage
bedtools multicov -bed bed/s_cerevisiae_genes.bed -bams bam/1M68_pH5_0.04CO2_R1.bam bam/1M69_pH5_0.04CO2_R2 > bedcoverage/s_cerevisiae_genes_multicov.tsv



bedtools multicov  -bed bed/s_cerevisiae_genes.bed -bams bam/1M68_pH5_0.04CO2_R1.bam bam/1M69_pH5_0.04CO2_R2.bam > bedcoverage/s_cerevisiae_genes_multicov.tsv