Lecture 7

Set up directory and copy over files

Check out the first few lines from the abridged results file

column -t fusion/kd3.fusion_predictions.abridged.tsv | less -S

Set up variables

module load StdEnv/2020 stringtie/2.1.5
awk '$1=="chr6"' gencode.gtf > gencode.chr6.gtf

---

Apply stringtie without and without reference annotations

stringtie --rf -o guided/ct1.gtf ct1.bam
stringtie --rf -G gencode.chr6.gtf -o denovo/ct1.gtf ct1.bam

---

Combine sample-specific outputs

stringtie --merge -o denovo.gtf denovo/*.gtf
stringtie --merge -o guided.gtf guided/*.gtf

---

Compare with reference annotations using gffcompare

module load nixpkgs/16.09 gffcompare/0.11.6
gffcompare -r gencode.chr6.gtf -o denovo denovo.gtf
gffcompare -r gencode.chr6.gtf -o guided guided.gtf

---

Check out the summary

less denovo.stats
less guided.stats

---

How many novel transcripts did we find?

awk '$3=="transcript"' guided.annotated.gtf| grep 'class_code "u"' | wc -l 

Set up environment

export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/r_env/bin:$PATH"

---

Parse out junctions with regtools

for b in *.bam; do regtools junctions extract -a 8 -m 50 -M 500000 -s RF $b -o $b.junc; done

---

Note sample information for downstream use

ls -A1 *junc > juncfiles
sed 's/.junc//' juncfiles \
  | awk -v OFS='\t' '{print $1, substr($1,0,2) == "ct" ? "control" : "knockdown"}' \
  > groupfiles

---

Run leafCutter

export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/r_env_V3/bin:$PATH"
python /home/hgen_share/utils/leafcutter/clustering/leafcutter_cluster_regtools.py -j juncfiles -m 50 -l 500000
Rscript /home/hgen_share/utils/leafcutter/scripts/leafcutter_ds.R \
  --num_threads 1 \
  --exon_file exons.txt.gz \
  --min_samples_per_intron 3 \
  leafcutter_perind_numers.counts.gz \
  groupfiles

---

Pack up results for visualization

Rscript /home/hgen_share/utils/leafcutter/leafviz/prepare_results.R \
  --meta_data_file groupfiles \
  --code KD \
  leafcutter_perind_numers.counts.gz \
  leafcutter_ds_cluster_significance.txt \
  leafcutter_ds_effect_sizes.txt \
  hg38 \
  -o viz.rda

---

Download viz.rda to your local computer

Set up directory and prepare files for input into rMATS

export PATH="/home/hgen_share/Anaconda/bin:$PATH"
export PATH="/home/hgen_share/Anaconda/envs/rmats/bin:$PATH"
mkdir output_rmats
mkdir tmp_rmats
echo "/home/hgen_share/lec7/ct1.bam,/home/hgen_share/lec7/ct2.bam,/home/hgen_share/lec7/ct3.bam" > control_list.txt
echo "/home/hgen_share/lec7/kd1.bam,/home/hgen_share/lec7/kd2.bam,/home/hgen_share/lec7/kd3.bam" > kd_list.txt

---

Run rMATS

python /home/hgen_share/Anaconda/envs/rmats/rMATS/rmats.py --b1 control_list.txt --b2 kd_list.txt --gtf /home/hgen_share/ref/hg38/ref-transcripts.gtf --od output_rmats -t paired --libType fr-firststrand --readLength 150 --nthread 1 --tmp tmp_rmats --variable-read-length

---