# Haemonchus contortus genome paper
## Section 5: Transcriptome complexity is generated by extensive cis- and trans- splicing

1. [Splice Leader analyses](#sl_analysis)
2. [LeafCutter](#leafcutter)


---
## Splice Leader analyses <a name="sl_analysis"></a>

```shell

cd /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME
mkdir SPLICE_LEADERS
cd SPLICE_LEADERS

ln -s ../../REF/HAEM_V4_final.chr.fa REF.fa
ln -s ../TRANSCRIPTOME_CURATION/HCON_V4_WBP11plus_190125.ips.gff3 ANNOTATION.gff3

```

Splice leader sequences
SL1 - GGTTTAATTACCCAAGTTTGAG
SL2 - GGTTTTAACCCAGTATCTCAAG


#### Run the splice leader analysis tools

```shell
# run splice leader finder scripts for SL1 and SL2 sequences
bsub.py --queue yesterday 1 merged_AdultF_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_AdultF_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultF_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultF_R2.fastq.gz
bsub.py --queue yesterday 1 merged_AdultM_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_AdultM_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultM_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultM_R2.fastq.gz
bsub.py --queue yesterday 1 merged_BXC_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_BXC_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_BXC_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_BXC_R2.fastq.gz
bsub.py --queue yesterday 1 merged_EXL3_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_EXL3_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_EXL3_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_EXL3_R2.fastq.gz
bsub.py --queue yesterday 1 merged_gut_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_gut_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_gut_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_gut_R2.fastq.gz
bsub.py --queue yesterday 1 merged_ISE_BXWF_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_ISE_BXWF_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_ISE_BXWF_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_ISE_BXWF_R2.fastq.gz
bsub.py --queue yesterday 1 merged_L1_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_L1_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L1_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L1_R2.fastq.gz
bsub.py --queue yesterday 1 merged_L4_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_L4_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L4_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L4_R2.fastq.gz
bsub.py --queue yesterday 1 merged_SHL3_SL1 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_SHL3_SL1 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTAATTACCCAAGTTTGAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_SHL3_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_SHL3_R2.fastq.gz


bsub.py --queue yesterday 1 merged_AdultF_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_AdultF_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultF_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultF_R2.fastq.gz
bsub.py --queue yesterday 1 merged_AdultM_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_AdultM_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultM_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_AdultM_R2.fastq.gz
bsub.py --queue yesterday 1 merged_BXC_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_BXC_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_BXC_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_BXC_R2.fastq.gz
bsub.py --queue yesterday 1 merged_EXL3_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_EXL3_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_EXL3_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_EXL3_R2.fastq.gz
bsub.py --queue yesterday 1 merged_gut_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_gut_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_gut_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_gut_R2.fastq.gz
bsub.py --queue yesterday 1 merged_ISE_BXWF_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_ISE_BXWF_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_ISE_BXWF_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_ISE_BXWF_R2.fastq.gz
bsub.py --queue yesterday 1 merged_L1_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_L1_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L1_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L1_R2.fastq.gz
bsub.py --queue yesterday 1 merged_L4_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_L4_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L4_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_L4_R2.fastq.gz
bsub.py --queue yesterday 1 merged_SHL3_SL2 ~sd21/bash_scripts/run_spliceleader_finder.sh merged_SHL3_SL2 $PWD/REF.fa $PWD/ANNOTATION.gff3 GGTTTTAACCCAGTATCTCAAG 10 /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_SHL3_R1.fastq.gz /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW/merged_SHL3_R2.fastq.gz
```
where "run_spliceleader_finder.sh" is:

```bash
#!/usr/bin/env bash
########################################################################

# SPLICE LEADER (SL) FINDER

########################################################################

# Usage: ~sd21/bash_scripts/run_spliceleader_finder.sh <PREFIX> <REFERENCE> <GFF> <splice leader sequence> <SL_MIN_LENGTH> <R1.fastq> <R2.fastq>

# Requirements in path
#--- samtools-1.6
#--- bedtools v2.17.0
#--- hisat2


# Notes
#--- use bsub to run
#--- relies on a Augustus-like GFF to detect overlap of SL with transcription start sites, and may fail to to the last step if not compatible. Can hack if needed.


# Modifications
#- 180712: added a minimum length of SL parameter for cutadapt to work with - originally set to 10 bp, but probably should be set to ~14 based on Alans script to find SLrelated seq by chance in genome
#- 180712: changed length of upstream window from start codon from 50 to 100 bp when looking for overlap between SL sequence and transcripts


########################################################################

PREFIX=$1
REFERENCE=$2
GFF=$3
SL_SEQ=$4
SL_MIN_LENGTH=$5
R1=$6
R2=$7


if [ "$#" -eq 0 ] || [ "$1" == "-h" ] || [ "$1" == "--help" ]; then
	echo ""
    echo "Usage: ~sd21/bash_scripts/run_spliceleader_finder.sh <PREFIX> <REFERENCE> <GFF> <splice leader sequence> <SL_MIN_LENGTH> <R1.fastq> <R2.fastq>"
    echo ""
    exit 0
fi

if [ -d "${PREFIX}_SL_ANALYSIS_out" ]; then
	echo -e "\nThere is already a directory started with this sample name. Rename and start again\n"
    exit 0
fi


mkdir ${PREFIX}_SL_ANALYSIS_out
cd ${PREFIX}_SL_ANALYSIS_out
ln -s ${R1}
ln -s ${R2}
ln -s ${REFERENCE}
ln -s ${GFF}


# cutadapt - http://cutadapt.readthedocs.io/en/stable/

# R1 reads with mates
cutadapt --trimmed-only --prefix=${PREFIX}_SL_trimmed_ --overlap=${SL_MIN_LENGTH} --error-rate=0.1 --info-file=info1.tmp -g ${SL_SEQ} -o ${PREFIX}.trimmed.1_1.tmp.fastq -p ${PREFIX}.trimmed.1_2.tmp.fastq ${R1} ${R2} > ${PREFIX}_cutadapt_trim.summary


# get the reverse reads by switching R1 and R2 around
cutadapt --trimmed-only --prefix=${PREFIX}_SL_trimmed_ --overlap=${SL_MIN_LENGTH} --error-rate=0.1 --info-file=info2.tmp -g ${SL_SEQ} -o ${PREFIX}.trimmed.2_2.tmp.fastq -p ${PREFIX}.trimmed.2_1.tmp.fastq ${R2} ${R1} >> ${PREFIX}_cutadapt_trim.summary


# extract adapter hit sequences, sort and count, make a sequence logo
#--- extract adapters
cat info1.tmp info2.tmp | awk '{if($2!="-1") print $0}' | cut -f 6 | sort | uniq -c | sort -k1,1n > ${PREFIX}_adapter_hits.counts
#--- make a fasta
while read count sequence; do echo -e ">adapter_${count}\n$sequence" >> ${PREFIX}_adapters.fa ; done < ${PREFIX}_adapter_hits.counts
#--- extend variable sequence to a fixed length (n=25bp) using "n"s to pad the sequence, and generate a sequence logo
awk '$1~">"{print $0}$1!~">"{tmp="";for(i=1;i<25-length($0)+1;i++){tmp=tmp"n"};print tmp""$0}' ${PREFIX}_adapters.fa >> ${PREFIX}_adapters.fixedlength.fa
~sd21/lustre118_link/software/bin/weblogo --format pdf --datatype fasta --ignore-lower-case < ${PREFIX}_adapters.fixedlength.fa > ${PREFIX}_adapters.fixedlength.pdf


# merge the two datasets together
cat ${PREFIX}.trimmed.1_1.tmp.fastq ${PREFIX}.trimmed.2_1.tmp.fastq > ${PREFIX}.trimmed.merged_R1.tmp.fastq
cat ${PREFIX}.trimmed.1_2.tmp.fastq ${PREFIX}.trimmed.2_2.tmp.fastq > ${PREFIX}.trimmed.merged_R2.tmp.fastq


# map reads using histat
/nfs/users/nfs_a/ar11/hisat2-2.0.0-beta/hisat2-build -p 1 ${REFERENCE} ${PREFIX}.REFidx.tmp
/nfs/users/nfs_a/ar11/hisat2-2.0.0-beta/hisat2 ${PREFIX}.REFidx.tmp -1 ${PREFIX}.trimmed.merged_R1.tmp.fastq -2 ${PREFIX}.trimmed.merged_R2.tmp.fastq -q -S ${PREFIX}.tmp.sam


# extract trimmed and mapped SL reads
samtools-1.6 view ${PREFIX}.tmp.sam -H > ${PREFIX}.tmp.sam.header
grep "${PREFIX}_SL_trimmed_" ${PREFIX}.tmp.sam > ${PREFIX}.tmp.sam.body
cat ${PREFIX}.tmp.sam.header ${PREFIX}.tmp.sam.body > ${PREFIX}.SL-only.tmp.sam
samtools-1.6 view -bS ${PREFIX}.SL-only.tmp.sam > ${PREFIX}.SL-only.tmp.bam
samtools-1.6 sort ${PREFIX}.SL-only.tmp.bam -o ${PREFIX}.SL-only.sorted.bam
samtools index -b ${PREFIX}.SL-only.sorted.bam


# extract coordinates, correcting for read direction, and output in bed format
bedtools bamtobed -i ${PREFIX}.SL-only.sorted.bam | grep "${PREFIX}_SL_trimmed_" | awk '{if($6=="+" && $5> 0) print $1,$2-1,$2+1,$4,$5,$6; else if($6=="-" && $5> 0) print $1,$3-1,$3+1,$4,$5,$6}' OFS="\t" > ${PREFIX}.SL-only.tmp.bed


# extract dinuclotide at splice site
bedtools getfasta -s -fi ${REFERENCE} -bed ${PREFIX}.SL-only.tmp.bed -fo ${PREFIX}.dinucleotide.tmp.fasta
cat ${PREFIX}.dinucleotide.tmp.fasta | grep -v ">" > ${PREFIX}.dinucleotide.tmp.seq


# extract twenty-mer at splice site
bedtools bamtobed -i ${PREFIX}.SL-only.sorted.bam | grep "${PREFIX}_SL_trimmed_" | awk '{if($6=="+" && $5> 0) print $1,$2-10,$2+10,$4,$5,$6; else if($6=="-" && $5> 0) print $1,$3-10,$3+10,$4,$5,$6}' OFS="\t" > ${PREFIX}.twentymer.tmp.bed
bedtools getfasta -s -fi  ${REFERENCE} -bed ${PREFIX}.twentymer.tmp.bed -fo ${PREFIX}.twentymer.fasta
cat ${PREFIX}.twentymer.fasta | grep -v ">" >  ${PREFIX}.twentymer.tmp.seq
~sd21/lustre118_link/software/bin/weblogo --format pdf < ${PREFIX}.twentymer.fasta > ${PREFIX}.twentymer.weblogo.pdf


# combine data
paste ${PREFIX}.SL-only.tmp.bed ${PREFIX}.dinucleotide.tmp.seq ${PREFIX}.twentymer.tmp.seq > ${PREFIX}.SL-only.bed


# calculate overlap with transcripts

awk '{if($3=="mRNA" || $3=="transcript") print $9}' ${GFF} | sed -e 's/ID=//g' -e 's/Parent=//g' -e 's/Name=//g' -e 's/;/\t/g' > ${PREFIX}.gene_transcript.list

#--- extract transcripts from GFF, and generate a flaking window around the start site 50 bp upstream and 20 downstream

while read -r transcript gene name other; do \
	grep "Parent=${transcript}" ${GFF} | \
	grep "CDS" | sort -k1,1 -k4,4n | \
	awk -v transcript_id="${transcript}" -v gene_id="${gene}" -v name_id="${name}" ' $7=="+" {print $1,$4-200,$4+30,"ID="transcript_id";Parent="gene_id";Name="name_id,".",$7; p=1 ; exit } END{if (p != 1) {print $1,$5-30,$5+200,"ID="transcript_id";Parent="gene_id";Name="name_id,".",$7}}' OFS="\t" >> ${PREFIX}_transcript_start_windows.tmp.bed ; \
	done < ${PREFIX}.gene_transcript.list

# clean up a little
#--- make sure all bed coords have stranded information - those that dont are from features without CDSs and therefore not needed
#--- fix and start coordinates in bed file that are less than 0. Replace with start coordinate of 1.

awk '{if($6=="+" || $6=="-") print $0}' ${PREFIX}_transcript_start_windows.tmp.bed | sed 's/-[0-9][0-9]*/1/' > ${PREFIX}_transcript_start_windows.bed

#--- calculate coverage of SL sequences in transcript start site windows
bedtools coverage -s -b ${PREFIX}_transcript_start_windows.bed -a ${PREFIX}.SL-only.bed | sort -k1,1 -k2,2n > ${PREFIX}_transcript_start_windows.SL.coverage

#--- extract internal CDS sequences from GFF, and generate a flaking window around the start site 10 bp upstream and 10 downstream
while read -r transcript gene name other; do \
	grep "Parent=${transcript}$" ${GFF} | \
	grep "CDS" | sort -k1,1 -k4,4n | \
	sed '1d;$d' | \
	awk -v transcript_id="${transcript}" -v gene_id="${gene}" -v name_id="${name}" '{print $1,$4-10,$5+10,"ID="transcript_id";Parent="gene_id";Name="name_id,".",$7}' OFS="\t" >> ${PREFIX}_internal_cds_windows.tmp.bed ; \
	done < ${PREFIX}.gene_transcript.list

awk '{if($6=="+" || $6=="-") print $0}' ${PREFIX}_internal_cds_windows.tmp.bed | sed 's/-[0-9][0-9]*/1/' > ${PREFIX}_internal_cds_windows.bed

bedtools coverage -s -b ${PREFIX}_internal_cds_windows.bed -a ${PREFIX}.SL-only.bed | sort -k1,1 -k2,2n > ${PREFIX}_internal_cds_windows.SL.coverage


#--- calculate coverage
awk '{if($7>1) print $4}' ${PREFIX}_transcript_start_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > ${PREFIX}_transcript_start_windows.SL.genelist

awk '{if($7>1) print $4}' ${PREFIX}_internal_cds_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 | sort | uniq -c > ${PREFIX}_internal_cds_windows.SL.genelist

# final cleanup
#rm *.tmp.*

```


```shell

awk '{if($6=="+" || $6=="-") print $0}' ${PREFIX}_transcript_start_windows.tmp.bed | sed 's/-[0-9][0-9]*/1/' transcript_start_windows.bed
```


### Run subsequence analyses to collate and summarise data
```shell
# collate all splice window coordinates
cat *SL1_SL_ANALYSIS_out/*SL-only.bed > SL1.SL-only.bed
cat *SL2_SL_ANALYSIS_out/*SL-only.bed > SL2.SL-only.bed

# get the transcript start and internal window coordinates - these are identical in each of the sub-analyses, so just choose one set
cp merged_EXL3_SL1_SL_ANALYSIS_out/merged_EXL3_SL1_transcript_start_windows.bed transcript_start_windows.bed
cp merged_EXL3_SL1_SL_ANALYSIS_out/merged_EXL3_SL1_internal_cds_windows.bed internal_cds_windows.bed

# get SL coverage per transcript start for both SL1 and SL2 sequences
bedtools coverage -s -b transcript_start_windows.bed -a SL1.SL-only.bed | sort -k1,1 -k2,2n > SL1_transcript_start_windows.SL.coverage
bedtools coverage -s -b transcript_start_windows.bed -a SL2.SL-only.bed | sort -k1,1 -k2,2n > SL2_transcript_start_windows.SL.coverage


# get SL coverage per internal CDSs (skipping terminal CDSs) for both SL1 and SL2 sequences
bedtools coverage -s -b internal_cds_windows.bed -a SL1.SL-only.bed | sort -k1,1 -k2,2n > SL1_internal_cds_windows.SL.coverage
bedtools coverage -s -b internal_cds_windows.bed -a SL2.SL-only.bed | sort -k1,1 -k2,2n > SL2_internal_cds_windows.SL.coverage


awk '{if($7>=1) print $4}' SL1_transcript_start_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > transcript_start_windows.SL1.genelist
awk '{if($7>=1) print $4}' SL1_transcript_start_windows.SL.coverage | wc -l
#>> 7895 transcripts
cut -f1 -d "-" transcript_start_windows.SL1.genelist | sort | uniq | wc -l
#>> 7293 genes



# awk '{if($7>=1) print $4}' SL1_internal_cds_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > internal_cds_windows.SL1.genelist
# awk '{if($7>=1) print $4}' SL1_internal_cds_windows.SL.coverage | wc -l
# #>> 5933 transcripts
# cut -f1 -d "-" internal_cds_windows.SL1.genelist | sort | uniq | wc -l
# #>> 3824 genes



awk '{if($7>=1) print $4}' SL2_transcript_start_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > transcript_start_windows.SL2.genelist
awk '{if($7>=1) print $4}' SL2_transcript_start_windows.SL.coverage | wc -l
#>> 1222 transcripts
cut -f1 -d "-" transcript_start_windows.SL2.genelist | sort | uniq | wc -l
#>> 1139 genes

# awk '{if($7>=1) print $4}' SL2_internal_cds_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > internal_cds_windows.SL2.genelist
# awk '{if($7>=1) print $4}' SL2_internal_cds_windows.SL.coverage | wc -l
# #>> 950 transcripts
# cut -f1 -d "-" internal_cds_windows.SL2.genelist | sort | uniq | wc -l
# #>> 766 genes


# #>> 3824 genes
```


### comparion SL counts per gene, to check coverage of genes that have both SL1 and SL2 sequences in reads
```bash
paste SL1_transcript_start_windows.SL.coverage SL2_transcript_start_windows.SL.coverage > test
awk '{if($7>0 || $17>0) print $7,$17}' OFS="\t" test > SL1v2_comparison.counts
```
```R
# load libraries
library(ggplot2)

# read data
a<-read.table("SL1v2_comparison.counts",header=F,sep="\t")

# make some plots
ggplot(a,aes(log10(V1),log10(V2)))+geom_hex(bins = 30)+labs(x="Number of SL1 reads [log10(count)]", y="Number of SL2 reads [log10(count)]")
ggplot(a,aes((V1+0.001),(V2+0.001)))+geom_hex(bins = 30)+labs(x="Number of SL1 reads [log10(count)]", y="Number of SL2 reads [log10(count)]")+scale_y_log10()+scale_x_log10()
```


### distance and density of genes, and their relationship between SL1 and SL2
```bash

awk -F '[\t=;]' '{if($3=="mRNA" && $7=="+") print $1,$4,$5,$10,$6,$7}' OFS="\t" ANNOTATION.gff3 | sort -k1,1 -k2,2n > transcripts.pos_strand.bed
awk -F '[\t=;]' '{if($3=="mRNA" && $7=="-") print $1,$4,$5,$10,$6,$7}' OFS="\t" ANNOTATION.gff3 | sort -k1,1 -k2,2n > transcripts.neg_strand.bed


bedtools-2 spacing -i transcripts.pos_strand.bed | sed 's/\.$/NA/g' > transcripts.pos_strand.spacing.bed
bedtools-2 spacing -i transcripts.neg_strand.bed | sed 's/\.$/NA/g' > transcripts.neg_strand.spacing.bed

awk '{print $4,$7}' OFS="\t" transcripts.pos_strand.spacing.bed > transcripts.pos_strand.spacing
awk '{print $4,$7}' OFS="\t" transcripts.neg_strand.spacing.bed > transcripts.neg_strand.spacing

# fix neg strand coords to correct for orientation
cut -f 1 transcripts.neg_strand.spacing | head -n -1 > neg.1.tmp
cut -f 2 transcripts.neg_strand.spacing | tail -n +2 > neg.2.tmp
paste neg.1.tmp neg.2.tmp > transcripts.neg_strand.spacing


awk '{print $1,"SL1"}' OFS="\t" transcript_start_windows.SL1.genelist > mRNA_SL1.list
awk '{print $1,"SL2"}' OFS="\t" transcript_start_windows.SL2.genelist > mRNA_SL2.list

cat mRNA_SL1.list mRNA_SL2.list > mRNA_SLall.list
```


### hybrid SL1/SL2 shared sites
```
bedtools-2 intersect -s -a SL1.SL-only.bed -b SL2.SL-only.bed > SL1SL2_sharedsites.bed

bedtools coverage -s -b transcript_start_windows.bed -a SL1SL2_sharedsites.bed | sort -k1,1 -k2,2n > SL1SL2_sharedsites_transcript_start_windows.SL.coverage

awk '{if($7>=1) print $4}' SL1SL2_sharedsites_transcript_start_windows.SL.coverage | sed -e 's/;/\t/g' -e 's/ID=//g' | cut -f1 > SL1SL2_sharedsites.genelist
awk '{if($7>=1) print $4}' SL1SL2_sharedsites_transcript_start_windows.SL.coverage | wc -l
# #>> 1021 transcripts
cut -f1 -d "-" SL1SL2_sharedsites.genelist | sort | uniq | wc -l
# #>> 1019 genes
```

### make some plots in R
```R
# load libraries
library(ggplot2)
library(patchwork)
library("ggsci")

# read data
pos<-read.table("transcripts.pos_strand.spacing",header=F)
neg<-read.table("transcripts.neg_strand.spacing",header=F)
sl<-read.table("mRNA_SLall.list",header=F)

# reformat data
data_pos<-dplyr::left_join(pos,sl,by="V1")
data_neg<-dplyr::left_join(neg,sl,by="V1")

sl1_pos<-data_pos[(data_pos$V2.y=="SL1"),]
sl2_pos<-data_pos[(data_pos$V2.y=="SL2"),]

sl1_neg<-data_neg[(data_neg$V2.y=="SL1"),]
sl2_neg<-data_neg[(data_neg$V2.y=="SL2"),]


par(mfrow=c(2,2))
plot(density(log10(sl1_pos$V2.x),na.rm=T),xlim=c(0,6))
plot(density(log10(sl1_neg$V2.x),na.rm=T),xlim=c(0,6))
plot(density(log10(sl2_pos$V2.x),na.rm=T),xlim=c(0,6))
plot(density(log10(sl2_neg$V2.x),na.rm=T),xlim=c(0,6))


plot_pos_all <- ggplot()+geom_histogram(aes(log10(data_pos$V2.x),fill=data_pos$V2.y),bins = 100, alpha = 0.5,position="identity")+theme_bw()+xlab("Distance to upstream gene (log10(bp))")+ylab("Count")+scale_fill_npg(na.value="grey90")+xlim(0,6)+ylim(0,500)
plot_neg_all <- ggplot()+geom_histogram(aes(log10(data_neg$V2.x),fill=data_neg$V2.y),bins = 100, alpha = 0.5,position="identity")+theme_bw()+xlab("Distance to upstream gene (log10(bp))")+ylab("Count")+scale_fill_npg(na.value="grey90")+xlim(0,6)+ylim(0,500)
plot_pos <- ggplot()+geom_histogram(aes(log10(sl1_pos$V2.x),stat(ndensity)),bins = 100,fill = "#E64B35B2", alpha = 0.5)+geom_histogram(aes(log10(sl2_pos$V2.x),stat(ndensity)),bins = 100,fill = "#4DBBD5B2", alpha = 0.5)+theme_bw()+xlab("Distance to upstream gene (log10(bp))")+ylab("Density")+xlim(0,6)
plot_neg <- ggplot()+geom_histogram(aes(log10(sl1_neg$V2.x),stat(ndensity)),bins = 100,fill = "#E64B35B2", alpha = 0.5)+geom_histogram(aes(log10(sl2_neg$V2.x),stat(ndensity)),bins = 100,fill = "#4DBBD5B2", alpha = 0.5)+theme_bw()+xlab("Distance to upstream gene (log10(bp))")+ylab("Density")+xlim(0,6)


# combined pos and neg
data_all <- dplyr::bind_rows(data_pos,data_neg)
sl1_all<-data_all[(data_all$V2.y=="SL1"),]
sl2_all<-data_all[(data_all$V2.y=="SL2"),]


p1 <- ggplot()+
     geom_histogram(aes(log10(data_all$V2.x),fill=(data_all$V2.y)),bins=100, alpha = 0.7,position="identity")+
     theme_bw()+
     xlab("Distance to upstream gene (log10(bp))")+ylab("Count")+scale_fill_npg(na.value="grey80")+
     xlim(0,6)+ylim(0,500)

p2 <- ggplot()+
     geom_histogram(aes(log10(sl1_all$V2.x),stat(ndensity)),bins = 100,fill = "#E64B35B2", alpha = 0.5)+
     geom_histogram(aes(log10(sl2_all$V2.x),stat(ndensity)),bins = 100,fill = "#4DBBD5B2", alpha=0.5)+
     theme_bw()+xlab("Distance to upstream gene (log10[bp])")+
     guides(fill = guide_legend(reverse = TRUE))+
     ylab("Density")+xlim(0,6)

p1 + p2 + plot_layout(ncol=2)
ggsave("SL_gene_distance.pdf")

```



### weblogo plots

```shell
*_adapters.fixedlength.fa
*_SL1.twentymer.fasta

cat *SL1_SL_ANALYSIS_out/*_adapters.fixedlength.fa > SL1_adapters_fixedlength.fa
cat *SL2_SL_ANALYSIS_out/*_adapters.fixedlength.fa > SL2_adapters_fixedlength.fa

~sd21/lustre118_link/software/bin/weblogo --format pdf --datatype fasta --ignore-lower-case <  SL1_adapters_fixedlength.fa >  SL1_adapters_fixedlength.weblogo.pdf
~sd21/lustre118_link/software/bin/weblogo --format pdf --datatype fasta --ignore-lower-case <  SL2_adapters_fixedlength.fa >  SL2_adapters_fixedlength.weblogo.pdf

cat *SL1_SL_ANALYSIS_out/*_SL1.twentymer.fasta > SL1_splicesite_20mer.fa
cat *SL2_SL_ANALYSIS_out/*_SL2.twentymer.fasta > SL2_splicesite_20mer.fa

~sd21/lustre118_link/software/bin/weblogo --format pdf --datatype fasta --ignore-lower-case <  SL1_splicesite_20mer.fa >  SL1_splicesite_20mer.weblogo.pdf
~sd21/lustre118_link/software/bin/weblogo --format pdf --datatype fasta --ignore-lower-case <  SL2_splicesite_20mer.fa >  SL2_splicesite_20mer.weblogo.pdf
```

[↥ **Back to top**](#top)

* * *


---
## Differential splicing w Leafcutter <a name="leafcutter"></a>

LEAFCUTTER
following http://davidaknowles.github.io/leafcutter/articles/Usage.html


NOTE - made modification to leafcutter scripts
- modified the following scripts
     - ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R
     - ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/prepare_results.R
     - ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/ds_plots.R
     - ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/run_leafviz.R

```shell
#changed the header from:
#  #!/usr/bin/env Rscript
#
# to
#
# #!/usr/bin/env /software/R-3.5.0/bin/Rscript
#
# to ensure the correct Rscript version was being used.
```

```shell
export PATH=$PATH:~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/

# raw data: /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/RAW

cd /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/LEAFCUTTER

ln -sf /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/TRANSCRIPTOME_CURATION/HCON_V4_WBP11plus_190125.ips.gff3 ANNOTATION.gff3
ln -s /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/REF/HAEM_V4_final.chr.fa REF.fa
```

Need to make a gtf of the annotation and curate it for leafcutter
```shell
# make gtf from gff
gffread -T -g REF.fa -o ANNOTATION.gtf ANNOTATION.gff3

# use leafcutter script to curate
perl ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/gtf2leafcutter.pl -o HCON_V4 ANNOTATION.gtf
```


```shell
# map RNAseq reads to reference using STAR
 while read name lane; do \
 		/nfs/users/nfs_s/sd21/lustre118_link/software/TRANSCRIPTOME/STAR/bin/Linux_x86_64/STAR \
 		--twopassMode Basic \
 		--runThreadN 8 \
 		--genomeDir /nfs/users/nfs_s/sd21/lustre118_link/REFERENCE_SEQUENCES/haemonchus_contortus/hc_v4chr_rnaseq_star_index \
 		--readFilesIn ../RAW/${lane}_1.fastq.gz ../RAW/${lane}_2.fastq.gz \
 		--readFilesCommand zcat \
 		--outSAMstrandField intronMotif \
 		--outSAMtype BAM Unsorted \
 		--outFileNamePrefix ${name}; \
 done < lanes.list


# make junc files
mkdir LEAFCUTTER_JUNCS

for bamfile in `ls STAR_MAPPING/*bam | awk -F "/" '{print $NF}' `; do \
		echo Converting $PWD/STAR_MAPPING/$bamfile to $PWD/LEAFCUTTER_JUNCS/$bamfile.junc; \
		sh bam2junc.sh $PWD/STAR_MAPPING/$bamfile $PWD/LEAFCUTTER_JUNCS/$bamfile.junc ; \
		echo $PWD/STAR_MAPPING/$bamfile.junc >> $PWD/LEAFCUTTER_JUNCS/test_juncfiles.txt; \
done



# make exon list
echo "chr start end strand gene_name" > HCON_V4.renamed.exons.leafcutter.list; \
     awk -F'[\t=;]' '{if($3=="exon") print $1,$4,$5,$7,$12}'  ANNOTATION.gff3 | sort -k1,1 -k 2,2n >> HCON_V4.renamed.exons.leafcutter.list
gzip HCON_V4.renamed.exons.leafcutter.list
```


### Run full analysis
```shell
# make comparison directories
mkdir LC_EGG_L1
mkdir LC_L1_SHL3
mkdir LC_EXL3_L4
mkdir LC_L4_AdultF
mkdir LC_EXL3_L4
mkdir LC_L4_AdultF
mkdir LC_L4_AdultM
mkdir LC_GUT_AdultF
mkdir LC_FEMALE_EGG
mkdir LC_SHL3_EXL3
mkdir LC_ADULTM_ADULTF


# copy junc files to each comparison dirctory
cp LEAFCUTTER_JUNCS/AdultF*junc LC_L4_AdultF/
cp LEAFCUTTER_JUNCS/AdultF*junc LC_GUT_AdultF
cp LEAFCUTTER_JUNCS/AdultF*junc LC_FEMALE_EGG
cp LEAFCUTTER_JUNCS/AdultM*junc LC_L4_AdultM
cp LEAFCUTTER_JUNCS/Adult*junc LC_ADULTM_ADULTF
cp LEAFCUTTER_JUNCS/egg*junc LC_EGG_L1
cp LEAFCUTTER_JUNCS/egg*junc LC_FEMALE_EGG
cp LEAFCUTTER_JUNCS/L1*junc LC_EGG_L1
cp LEAFCUTTER_JUNCS/L1*junc LC_L1_SHL3
cp LEAFCUTTER_JUNCS/L4*junc LC_EXL3_L4
cp LEAFCUTTER_JUNCS/L4*junc LC_L4_AdultF
cp LEAFCUTTER_JUNCS/L4*junc LC_L4_AdultM/
cp LEAFCUTTER_JUNCS/gut*junc LC_GUT_AdultF

cp LEAFCUTTER_JUNCS/SHL3*junc LC_SHL3_EXL3
cp LEAFCUTTER_JUNCS/EXL3*junc LC_SHL3_EXL3

# fix for renaming of SHL3 and EXL3 as per kallisto analysis
cd LC_SHL3_EXL3
mv SHL3_3_236476_2305_3914314Aligned.out.bam.junc EXL3_kallisto_7062_6_12_out.junc
mv EXL3_2_236476_2305_3914310Aligned.out.bam.junc SHL3_kallisto_7062_6_8_out.junc
rm EXL3_1_236476_2305_3914309Aligned.out.bam.junc

# move correct samples to other directories
cp EXL3*junc ../LC_EXL3_L4
cp SHL3*junc ../LC_L1_SHL3

cd ../


# make list of junc files in each comparison directory
for i in LC_* ; do \
		ls -1  $i/*junc | \
		xargs -n1 basename > ${i}/juncfiles.list; \
		done


# cluster splice junctions
#--- min depth 5
for i in LC_* ; do \
     cd ${i} && \
     python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 5 -o ${i}_cov5 -l 500000 && \
     cd ../; \
     done

#--- min depth 10
for i in LC_* ; do \
     cd ${i} && \
     python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 10 -o ${i}_cov10 -l 500000 && \
     cd ../; \
     done

#--- min depth 30 (default)
for i in LC_* ; do \
          cd ${i} && \
          python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 30 -o ${i}_cov30 -l 500000 && \
          cd ../; \
          done



# all samples
mkdir LEAFCUTTER_JUNCS_V4
cp LC_*/*junc LEAFCUTTER_JUNCS_V4
cd LEAFCUTTER_JUNCS_V4

ls -1 *junc > juncfiles.list

python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 5 -o all_samples_cov5 -l 500000

python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 10 -o all_samples_cov10 -l 500000

python2.7 ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/clustering/leafcutter_cluster.py -j juncfiles.list -m 30 -o all_samples_cov30 -l 500000


cd ..
```


### perform differential analysis
```
#--- made "samples_groups.list" files for each comparison manually, as the sample names were not really parsable.
#--- min depth 5
for i in LC_FEMALE_EGG LC_GUT_AdultF LC_L1_SHL3 LC_L4_AdultF LC_L4_AdultM LC_ADULTM_ADULTF LC_EGG_L1; do \
     		cd ${i} && \
     		~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
     		--min_samples_per_intron=3 --min_samples_per_group=3 --min_coverage=10 --exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov5_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov5 && \
     		cd ../ ; \
     		done

for i in LC_SHL3_EXL3 LC_EXL3_L4 ; do \
               cd ${i} && \
               ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
               --min_samples_per_intron=2 --min_samples_per_group=2 --min_coverage=10 --exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov5_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov5 && \
               cd ../ ; \
               done



#--- min depth 10
for i in LC_FEMALE_EGG LC_GUT_AdultF LC_L1_SHL3 LC_L4_AdultF LC_L4_AdultM LC_ADULTM_ADULTF LC_EGG_L1; do \
          cd ${i} && \
          ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
          --min_samples_per_intron=3 --min_samples_per_group=3 --min_coverage=10 --exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov10_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov10 && \
          cd ../ ; \
          done

for i in LC_SHL3_EXL3 LC_EXL3_L4 ; do \
          cd ${i} && \
          ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
          --min_samples_per_intron=2 --min_samples_per_group=2 -exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov10_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov10 && \
          cd ../ ; \
          done


#--- min depth 30
for i in LC_FEMALE_EGG LC_GUT_AdultF LC_L1_SHL3 LC_L4_AdultF LC_L4_AdultM LC_ADULTM_ADULTF LC_EGG_L1; do \
                    cd ${i} && \
                    ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
                    --min_samples_per_intron=3 --min_samples_per_group=3 --exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov30_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov30 && \
                    cd ../ ; \
                    done

for i in LC_SHL3_EXL3 LC_EXL3_L4 ; do \
                    cd ${i} && \
                    ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/scripts/leafcutter_ds.R \
                    --min_samples_per_intron=2 --min_samples_per_group=2 --exon_file=../HCON_V4_all_exons.txt.gz ${i}_cov30_perind_numers.counts.gz samples_groups.list --output_prefix=${i}_cov30 && \
                    cd ../ ; \
                    done


#--- made "samples_groups.list" files for each comparison manually, as the sample names were not really parsable
for i in LC*; do cd $i; ls -1 *junc | awk -F'[_]' '{print $0,$1}' OFS="\t" > samples_groups.list; cd ../ ;  done


# prepare results for visualisation
#--- min depth 5
for i in LC_* ; do \
		cd ${i} && \
		~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/prepare_results.R ${i}_cov5_perind_numers.counts.gz ${i}_cov5_cluster_significance.txt ${i}_cov5_effect_sizes.txt ../HCON_V4 --meta_data_file samples_groups.list --output=${i}_leafviz_min5.RData && \
		cd ../ ; \
		done

#--- min depth 10
for i in LC_* ; do \
		cd ${i} && \
		~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/prepare_results.R ${i}_cov10_perind_numers.counts.gz ${i}_cov5_cluster_significance.txt ${i}_cov5_effect_sizes.txt ../HCON_V4 --meta_data_file samples_groups.list --output=${i}_leafviz_min10.RData && \
		cd ../ ; \
		done



#--- min depth 30
for i in LC_* ; do \
          cd ${i} && \
          ~sd21/lustre118_link/software/TRANSCRIPTOME/leafcutter/leafviz/prepare_results.R ${i}_cov30_perind_numers.counts.gz ${i}_cov30_cluster_significance.txt ${i}_cov30_effect_sizes.txt ../HCON_V4 --meta_data_file samples_groups.list --output=${i}_leafviz_min30.RData && \
          cd ../ ; \
          done


# copy data to laptop for visualisation

for i in LC_* ; do cp ${i}/${i}_leafviz_min30.RData ~sd21/desktop_link/${i}_leafviz_min30.RData; done


#- on laptop, move data here: /Users/sd21/Documents/Work/aaa_projects/haemonchus_contortus/project_hc_genome/data/Section_4/leafcutter_data
cp /Volumes/sd21/desktop_link/*.RData .

# obviosuly needs leafcutter installed locally, which is ehre: http://davidaknowles.github.io/leafcutter/articles/Installation.html

# also copied

# to visualise on laptop
cd /Users/sd21/Documents/Work/aaa_projects/haemonchus_contortus/project_hc_genome/data/Section_4/leafcutter_data
cp ~/Documents/bioinformatics_bin/leafcutter/leafviz/server.R .
cp ~/Documents/bioinformatics_bin/leafcutter/leafviz/ui.R .

~/Documents/bioinformatics_bin/leafcutter/leafviz/run_leafviz.R LF_L4_AdultM_leafviz.RData
```



### Meta analysis - compare differential intro usage across all replicates for all lifestages
```bash
# working dir:
cd /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/LEAFCUTTER
mkdir METAANALYSIS_V4

# get sig clusters from  cluster_significance.txt  - top 50 based of FDR pvalue
for i in LC_* ; do \
          cd ${i} && \
          awk -F'[:\t]' '{if($7<0.05) print $2,$7}' OFS="\t" ${i}_cov30_cluster_significance.txt | sort -k2 -g | head -n 100 > ../METAANALYSIS_V4/${i}_sigclusters.list && \
          cd ../ ; \
          done


# extract introns from clusters that are signficant from effect_sizes.txt
rm METAANALYSIS_V4/sigclusters_introncoords.list

for i in LC_* ; do \
          cd ${i} && \
          while read cluster sig; do \
          grep $cluster ${i}_cov5_effect_sizes.txt; done < ../METAANALYSIS_V4/${i}_sigclusters.list  | awk -F'[:\t]' '{if($4>=5 || $4<=-5) print $1,$2,$3}' OFS=":" >> ../METAANALYSIS_V4/sigclusters_introncoords.list && \
          cd ../ ; \
          done

cd METAANALYSIS_V4
cat sigclusters_introncoords.list | sort -k1n | uniq > sigclusters_introncoords.sorted.list


# extract counts using cluster, intro coords
while read coords; do grep $coords <(zcat ../LEAFCUTTER_JUNCS_V4/all_samples_cov30_perind.counts.gz) ; done < sigclusters_introncoords.sorted.list > sigclusters_introncoords.junccounts


# ---- Sample columns in all_samples_cov30_perind.counts.gz
eggs1_236476_3881079Aligned.out.bam 20
eggs2_236476_3881080Aligned.out.bam 19
eggs3_236476_3881081Aligned.out.bam 14

L1_1_236476_3881082Aligned.out.bam 24
L1_2_236476_3881083Aligned.out.bam 10
L1_3_236476_3881084Aligned.out.bam 4

SHL3_1_236476_2305_3914312Aligned.out.bam 6
SHL3_2_236476_2305_3914313Aligned.out.bam 12
SHL3_kallisto_7062_6_8_out 22

EXL3_3_236476_2305_3914311Aligned.out.bam 5
EXL3_kallisto_7062_6_12_out 18

L4_1_236476_3881085Aligned.out.bam 17
L4_2_236476_3881086Aligned.out.bam 21
L4_3_236476_3881087Aligned.out.bam 11


AdultF1_236476_3881088Aligned.out.bam 2
AdultF2_236476_3881089Aligned.out.bam 23
AdultF3_236476_3881090Aligned.out.bam 8

AdultM1_236476_2305_3914303Aligned.out.bam 7
AdultM2_236476_2305_3914304Aligned.out.bam 9
AdultM3_236476_2305_3914305Aligned.out.bam 3

gut3_236476_635J_3914317Aligned.out.bam 13
gut2_236476_1589_3914316Aligned.out.bam 15
gut1_236476_1517_3914315Aligned.out.bam 16



# reorder the columns to fit the haem lifecycle, and then calc the splice frequency
awk '{print $20,$19,$14,$24,$10,$4,$6,$12,$22,$5,$18,$17,$21,$11,$2,$23,$8,$7,$9,$3,$13,$15,$16}' sigclusters_introncoords.junccounts |\
	sed 's/\//\t/g' |\
	awk '{print $1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001),$13/($14+=0.0001),$15/($16+=0.0001),$17/($18+=0.0001),$19/($20+=0.0001),$21/($22+=0.0001),$23/($24+=0.0001),$25/($26+=0.0001),$27/($28+=0.0001),$29/($30+=0.0001),$31/($32+=0.0001),$33/($34+=0.0001),$35/($36+=0.0001),$37/($38+=0.0001),$39/($40+=0.0001),$41/($42+=0.0001),$43/($44+=0.0001),$45/($46+=0.0001)}' OFS="\t" > sigclusters_introncoords.juncfreq
	awk '{if(($1+$2+$3+$4+$5+$6+$7+$8+$9+$10+$11+$12+$13+$14+$15+$16+$17+$18+$19+$20+$21+$22+$23)>=0.5) print $0}' OFS="\t" sigclusters_introncoords.juncfreq > sigclusters_introncoords.juncfreq2
```

# R commands to make heatmap
```R
#load libraries
library(gplots)
library(RColorBrewer)

pdf("allsamples.leafcutter_ds.top100.pdf")
a<-read.table("sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("Egg_1","Egg_2","Egg_3","L1_1","L1_2","L1_3","SHL3_1","SHL3_2","SHL3_3","EXL3_1","EXL3_2","L4_1","L4_2","L4_3","Adult_Female_1","Adult_Female_2","Adult_Female_3","Adult_Male_1","Adult_Male_2","Adult_Male_3","Gut_1","Gut_2","Gut_3")
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```


### repeat extract significant clusters for each pairwise comparison
```bash
# working dir
cd /nfs/users/nfs_s/sd21/lustre118_link/hc/GENOME/TRANSCRIPTOME/LEAFCUTTER


# get sig clusters from  cluster_significance.txt  - top 100 based of FDR pvalue
# --  done in the meta analysis step above

# extract intron coords per sample, sort
for i in LC_* ; do \
		cd ${i} && \
		awk -F'[:\t]' '{if($6<0.05) print $2,$7,$8}' OFS="\t" ${i}_cov30_cluster_significance.txt | sort -k2 -g | head -n 100 > ${i}_sigclusters.list
		while read cluster sig mrna_id; do grep $cluster ${i}_cov30_effect_sizes.txt; done < ${i}_sigclusters.list  | awk -F'[:\t]' '{if($4>=5 || $4<=-5) print $1,$2,$3}' OFS=":" >  ${i}_sigclusters_introncoords.list
		cat ${i}_sigclusters_introncoords.list | sort -k1n | uniq > ${i}_sigclusters_introncoords.sorted.list
		while read coords; do grep $coords <(zcat ${i}_cov30_perind.counts.gz) ; done < ${i}_sigclusters_introncoords.sorted.list > ${i}_sigclusters_introncoords.junccounts
		cd ../ ; \
		done
```

### compare pairs of lifestages for differential splicing

```bash
# eggs to L1
cd ../LC_EGG_L1
awk '{print $6,$5,$4,$3,$2,$7,$1}' LC_EGG_L1_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_EGG_L1_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_EGG_L1_sigclusters_introncoords.juncfreq; done < LC_EGG_L1_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_EGG_L1_sigclusters_introncoords.juncfreq > LC_EGG_L1_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_EGG_L1.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_EGG_L1_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("Egg_1","Egg_2","Egg_3","L1_1","L1_2","L1_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```

```bash
# L1 to SHL3
cd ../LC_L1_SHL3
awk '{print $5,$2,$3,$6,$4,$7,$1}' LC_L1_SHL3_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_L1_SHL3_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_L1_SHL3_sigclusters_introncoords.juncfreq ; done < LC_L1_SHL3_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_L1_SHL3_sigclusters_introncoords.juncfreq > LC_L1_SHL3_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_L1_SHL3.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_L1_SHL3_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("L1_1","L1_2","L1_3","SHL3_1","SHL3_2","SHL3_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```



```bash
# SHL3 to EXL3
cd ../LC_SHL3_EXL3
awk '{print $2,$3,$6,$4,$5,$1}' LC_SHL3_EXL3_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $14,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001)}' OFS="\t" > LC_SHL3_EXL3_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_SHL3_EXL3_sigclusters_introncoords.juncfreq ; done < LC_SHL3_EXL3_sigclusters.list
awk '{if(($2+$3+$4+$5+$6)>=0.5) print $2,$3,$4,$5,$6}' OFS="\t" LC_SHL3_EXL3_sigclusters_introncoords.juncfreq > LC_SHL3_EXL3_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_SHL3_EXL3.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_SHL3_EXL3_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("SHL3_1","SHL3_2","SHL3_3","EXL3_1","EXL3_2")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```


```bash
# EXL3 to L4
cd ../LC_EXL3_L4
awk '{print $4,$5,$2,$6,$3,$1}' LC_EXL3_L4_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $14,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001)}' OFS="\t" > LC_EXL3_L4_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_EXL3_L4_sigclusters_introncoords.juncfreq ; done < LC_EXL3_L4_sigclusters.list
awk '{if(($2+$3+$4+$5+$6)>=0.5) print $2,$3,$4,$5,$6}' OFS="\t" LC_EXL3_L4_sigclusters_introncoords.juncfreq > LC_EXL3_L4_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_EXL3_L4.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_EXL3_L4_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("EXL3_1","EXL3_2","L4_1","L4_2","L4_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```



```bash
# L4 to AdultMale
cd ../LC_L4_AdultM
awk '{print $2,$5,$3,$6,$7,$4,$1}' LC_L4_AdultM_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_L4_AdultM_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_L4_AdultM_sigclusters_introncoords.juncfreq ; done < LC_L4_AdultM_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_L4_AdultM_sigclusters_introncoords.juncfreq > LC_L4_AdultM_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_L4_AdultM.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_L4_AdultM_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("L4_1","L4_2","L4_3","Adult_Male_1","Adult_Male_2","Adult_Male_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```
cd ../



```bash
# L4 to AdultFemale
cd ../LC_L4_AdultF
awk '{print $2,$6,$4,$3,$5,$7,$1}' LC_L4_AdultF_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_L4_AdultF_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_L4_AdultF_sigclusters_introncoords.juncfreq ; done < LC_L4_AdultF_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_L4_AdultF_sigclusters_introncoords.juncfreq > LC_L4_AdultF_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_L4_AdultF.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_L4_AdultF_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("L4_1","L4_2","L4_3","Adult_Female_1","Adult_Female_2","Adult_Female_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))		
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```


```bash
# AdultMale to AdultFemale
cd ../LC_ADULTM_ADULTF
awk '{print $5,$7,$3,$2,$4,$6,$1}' LC_ADULTM_ADULTF_sigclusters_introncoords.junccounts |\
     sed -e 's/\//\t/g' -e 's/:/\t/g' |\
     awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq ; done < LC_ADULTM_ADULTF_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq > LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq2

while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}_${cluster}/g" LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq ; done < LC_ADULTM_ADULTF_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print NR"_"$1,$2,$3,$4,$5,$6,$7}' OFS="\t" LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq > LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_ADULTM_ADULTF.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_ADULTM_ADULTF_sigclusters_introncoords.juncfreq2",header=F,row.names=1)
colnames(a)<- c("Adult_Male_1","Adult_Male_2","Adult_Male_3","Adult_Female_1","Adult_Female_2","Adult_Female_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```


```bash
# AdultFemale to eggs
cd ../LC_FEMALE_EGG
awk '{print $2,$3,$7,$6,$5,$4,$1}' LC_FEMALE_EGG_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_FEMALE_EGG_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_FEMALE_EGG_sigclusters_introncoords.juncfreq ; done < LC_FEMALE_EGG_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_FEMALE_EGG_sigclusters_introncoords.juncfreq > LC_FEMALE_EGG_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_FEMALE_EGG.leafcutter_ds.top100.pdf")
library(gplots)
library(RColorBrewer)
a<-read.table("LC_FEMALE_EGG_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("Adult_Female_1","Adult_Female_2","Adult_Female_3","Egg_1","Egg_2","Egg_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```

cd ../

```bash
# AdultFemale to gut
cd ../LC_GUT_AdultF
awk '{print $2,$3,$7,$6,$5,$4,$1}' LC_GUT_AdultF_sigclusters_introncoords.junccounts |\
sed -e 's/\//\t/g' -e 's/:/\t/g' |\
awk '{print $16,$1/($2+=0.0001),$3/($4+=0.0001),$5/($6+=0.0001),$7/($8+=0.0001),$9/($10+=0.0001),$11/($12+=0.0001)}' OFS="\t" > LC_GUT_AdultF_sigclusters_introncoords.juncfreq
while read cluster sig mrna; do sed -i "s/${cluster}/${mrna}/g" LC_GUT_AdultF_sigclusters_introncoords.juncfreq ; done < LC_GUT_AdultF_sigclusters.list
awk '{if(($2+$3+$4+$5+$6+$7)>=0.5) print $2,$3,$4,$5,$6,$7}' OFS="\t" LC_GUT_AdultF_sigclusters_introncoords.juncfreq > LC_GUT_AdultF_sigclusters_introncoords.juncfreq2
```

```R-3.5.0
pdf("LC_GUT_AdultF.leafcutter_ds.top100.pdf",)
library(gplots)
library(RColorBrewer)
a<-read.table("LC_GUT_AdultF_sigclusters_introncoords.juncfreq2",header=F)
colnames(a)<- c("Adult_Female_1","Adult_Female_2","Adult_Female_3","Gut_1","Gut_2","Gut_3")
#heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,margins = c(10, 10))
heatmap.2(as.matrix(a),Colv=NULL,trace="none",col= colorRampPalette(brewer.pal(8, "Blues"))(25),labRow=FALSE,labCol=FALSE,margins = c(10, 10),dendrogram='none',density.info='none')
dev.off()
quit(save = "no")
```



### calculated summary stats
``` bash
# calculated summary stats
cd ../
for i in LC_*; do echo ${i}; \
# total genes
echo "total genes"; \
cut -f7 ${i}/*cov30_cluster_significance.txt | sort | uniq | wc -l ; \
echo "total significant genes"; \
awk '{if($6<0.05) print $7}' ${i}/*cov30_cluster_significance.txt | sort | uniq | wc -l ; \
# total introns
echo "total introns"; \
cat ${i}/*cov30_cluster_significance.txt | wc -l ; \
# total significant introns
echo "total significant introns"; \
awk '{if($6<0.05) print $0}' ${i}/*cov30_cluster_significance.txt | wc -l; \
echo " "
done


# calculate the number of genes in total impacted by differential splicing
cat LC_*/*cov30_cluster_significance.txt | awk '{if($6<0.05) print }' | cut -f7 | sed 's/,/\n/g' | sort | uniq | wc -l
#>1055
```
