Source code to produce the PSI plots

<h1 id="source-code-to-produce-the-psi-plots">Source code to produce the PSI plots</h1>
<p>The scripts were used to generate the PSI plots to appear in Radke, M and Badillo-Lisakowski, V et al &quot;Therapeutic inhibition of RBM20 improves diastolic function in a murine heart failure model and human engineered heart tissue&quot;.
deposited for reproducibility purposes.</p>
<p>This repository was updated on 27.01.2022 to include the <code>Snakefile</code> and <code>compute_PSI.pl</code>.</p>
<p>The workflow is based on Snakemake, see <code>Snakefile</code> for more information.
Workflow was called with:</p>
<pre><code class="lang-bash">snakemake
  -<span class="ruby">-configfile config.yml
</span>  -<span class="ruby">-printshellcmds
</span>  -<span class="ruby">-cluster <span class="hljs-string">"sbatch
</span></span>  -<span class="ruby"><span class="hljs-string">-mem {cluster.mem}
</span></span>  -<span class="ruby"><span class="hljs-string">-out {cluster.out}
</span></span>  -<span class="ruby"><span class="hljs-string">-error {cluster.out}
</span></span>  -<span class="ruby"><span class="hljs-string">c {cluster.c}
</span></span>  -<span class="ruby"><span class="hljs-string">-parsable "</span>
</span>  -<span class="ruby">-cluster-config cluster.yaml
</span>  -<span class="ruby">-jobs <span class="hljs-number">10</span></span>
</code></pre>
<h2 id="dependencies-">Dependencies:</h2>
<ul>
<li>stringtie/1.3.5</li>
<li>cufflinks/2.2.1</li>
<li>python2</li>
<li>R/3.5.1</li>
<li>perl</li>
<li>bedtools</li>
</ul>
<p>Dependencies were handled with Environment Modules on the Snakefile.</p>
<h2 id="workflow-steps-">Workflow steps:</h2>
<ol>
<li><p>Guided transcriptome assembly with StringTie:</p>
</li>
<li><p>Parameters:</p>
<ul>
<li>min_isoform_fraction: 0.2</li>
<li>min_anchor_length: 8</li>
<li>min_junction_coverage: 3</li>
<li>reference: /biodb/genomes/mus_musculus/GRCm38_85/GRCm38.85.gtf
The modified GRCm38 84 annotation accopanies this repository. The version of Ensembl has been used to study Titin and match the exons to domains.</li>
</ul>
</li>
<li><p>Merge annotation with cuffmerge:</p>
</li>
<li><p>Parameters:</p>
<ul>
<li>reference: /biodb/genomes/mus_musculus/GRCm38_85/GRCm38.85.gtf</li>
<li>min_isoform_fraction: 0.2</li>
</ul>
</li>
<li><p>Find &#39;containers&#39;, loci with that contains alternative splicing events, and calculate PSI from STAR <code>SJ.out.tab</code>  intron counts with compute_PSI.pl
PSI is define by the  ratio between the number of reads including an exon and total reads for that event (reads including + reads excluding).</p>
</li>
<li><p>Produce the plots using PSI_plot_aso.R script:
One plot per gene, only for gene with two or more exons, in addiiton:</p>
</li>
<li>Removde duplicated exons</li>
<li>Merged overlapping exons</li>
</ol>
<p>The accompanying test data reproduces figure 4.b of the paper:
<code>Rscript testdata/PSI_plot_aso.R testdata/output/ testdata/ldb3.gtf testdata/</code></p>
<p>The subset of the data was prepared with:</p>
<pre><code class="lang-{fish}"><span class="hljs-built_in">set</span> gtf <span class="hljs-string">"/beegfs/prj/Mitch_Gotthardt/Thiago/psi_plots/Mus_musculus.GRCm38.84_cleaned.gtf"</span>
grep <span class="hljs-string">'gene_name "Ldb3";'</span> <span class="hljs-variable">$gtf</span> | sort -k1,1 -k4,4n | bedtools merge <span class="hljs-_">-d</span> 1000000 <span class="hljs-_">-s</span> -c 5,6,7 -o first &gt; testdata/ldb3.bed
intersectBed -b testdata/ldb3.bed <span class="hljs-_">-a</span> <span class="hljs-variable">$gtf</span> -wa <span class="hljs-_">-s</span> &gt; testdata/ldb3.gtf

<span class="hljs-keyword">for</span> f <span class="hljs-keyword">in</span> stringtie_guided/*/psi.txt;
  mkdir testdata/(basename (dirname <span class="hljs-variable">$f</span>));
  intersectBed -b testdata/ldb3.bed <span class="hljs-_">-a</span> <span class="hljs-variable">$f</span> -wa <span class="hljs-_">-s</span> &gt; testdata/(basename (dirname <span class="hljs-variable">$f</span>))/psi.txt;
end
</code></pre>