

In [1]:

    
ls -lah ../data/









    



total 191M
drwxrwxr-x 2 ilya ilya 4.0K Sep 30 18:05 ./
drwxrwxr-x 5 ilya ilya 4.0K Sep 16 13:23 ../
-rw-rw-r-- 1 ilya ilya  24M Sep 30 08:54 7E2.pileup.gz
-rw-rw-r-- 1 ilya ilya  21M Sep 30 15:49 AAGCTA_R2.pileup.gz
-rw-rw-r-- 1 ilya ilya 146M Sep 29 12:04 BJ-HSR1_R1.fastq.gz
-rw-rw-r-- 1 ilya ilya 500K Sep 16 09:03 contigs.fasta
-rw-r----- 1 ilya ilya  644 Sep 16 09:00 dHSR1.fa
-rw-rw-r-- 1 ilya ilya 5.4K Sep 23 15:58 gradtimes.txt
-rw-rw-r-- 1 ilya ilya  445 Sep 16 09:00 hHSR-435.fa
-rw-rw-r-- 1 ilya ilya  611 Sep 16 09:00 hHSR.fa
-rw-rw-r-- 1 ilya ilya  126 Sep 16 09:04 rose.fa



In [2]:

    
!head ../data/rose.fa









    



>ROSE1
GCCGCCGAGAGGUGGCGUCCCCGACGCCUCAUGGGUCGGGAACGACUGAGACGGGCACCG
GUCGUGUCCGGUGCCGCUCGUAUCCAUUUUGCUCCUUGGAGGAUUUGGCUAUGCGCA



In [3]:

    
%%bash

cd ../data/
RNAfold -p -d2 --noPS --noLP -T 37 < rose.fa
cd -









    



>ROSE1
GCCGCCGAGAGGUGGCGUCCCCGACGCCUCAUGGGUCGGGAACGACUGAGACGGGCACCGGUCGUGUCCGGUGCCGCUCGUAUCCAUUUUGCUCCUUGGAGGAUUUGGCUAUGCGCA
((((((....))))))(((((((((.((....))))))))...))).......((((((((......))))))))((.((((.(((..((.(((....)))))..))).)))).)). (-57.30)
((((((....))))))(({((((((.((....)))))))).,,}}}.......((((((((......))))))))((.((((.(((..({.(((....)))})..))).)))).)). [-58.22]
((((((....))))))(((((((((.((....))))))))...))).......((((((((......))))))))((.((((.(((..((.(((....)))))..))).)))).)). {-57.30 d=5.20}
 frequency of mfe structure in ensemble 0.225946; ensemble diversity 7.63  
/home/ilya/src/biodata/sessions



In [4]:

    
ls -lah ../data/









    



total 191M
drwxrwxr-x 2 ilya ilya 4.0K Oct 14 09:18 ./
drwxrwxr-x 5 ilya ilya 4.0K Sep 16 13:23 ../
-rw-rw-r-- 1 ilya ilya  24M Sep 30 08:54 7E2.pileup.gz
-rw-rw-r-- 1 ilya ilya  21M Sep 30 15:49 AAGCTA_R2.pileup.gz
-rw-rw-r-- 1 ilya ilya 146M Sep 29 12:04 BJ-HSR1_R1.fastq.gz
-rw-rw-r-- 1 ilya ilya 500K Sep 16 09:03 contigs.fasta
-rw-r----- 1 ilya ilya  644 Sep 16 09:00 dHSR1.fa
-rw-rw-r-- 1 ilya ilya 5.4K Sep 23 15:58 gradtimes.txt
-rw-rw-r-- 1 ilya ilya  445 Sep 16 09:00 hHSR-435.fa
-rw-rw-r-- 1 ilya ilya  611 Sep 16 09:00 hHSR.fa
-rw-rw-r-- 1 ilya ilya  11K Oct 14 09:18 ROSE1_dp.ps
-rw-rw-r-- 1 ilya ilya  126 Sep 16 09:04 rose.fa



In [5]:

    
!head -n 25 ../data/ROSE1_dp.ps









    



%!PS-Adobe-3.0 EPSF-3.0
%%Title: RNA Dot Plot
%%Creator: PS_dot.c,v 1.38 2007/02/02 15:18:13 ivo Exp $, ViennaRNA-2.1.9
%%CreationDate: Wed Oct 14 09:18:15 2015
%%BoundingBox: 66 211 518 662
%%DocumentFonts: Helvetica
%%Pages: 1
%%EndComments

%Options: -noLP -d2 
% 
%This file contains the square roots of the base pair probabilities in the form
% i  j  sqrt(p(i,j)) ubox

%%BeginProlog
/DPdict 100 dict def
DPdict begin
/logscale false def
/lpmin 1e-05 log def

/box { %size x y box - draws box centered on x,y
   2 index 0.5 mul sub            % x -= 0.5
   exch 2 index 0.5 mul sub exch  % y -= 0.5
   3 -1 roll dup rectfill
} bind def



In [6]:

    
!tail -n 25 ../data/ROSE1_dp.ps









    



54 75 0.9500000 lbox
55 74 0.9500000 lbox
56 73 0.9500000 lbox
57 72 0.9500000 lbox
58 71 0.9500000 lbox
59 70 0.9500000 lbox
60 69 0.9500000 lbox
61 68 0.9500000 lbox
76 116 0.9500000 lbox
77 115 0.9500000 lbox
79 113 0.9500000 lbox
80 112 0.9500000 lbox
81 111 0.9500000 lbox
82 110 0.9500000 lbox
84 108 0.9500000 lbox
85 107 0.9500000 lbox
86 106 0.9500000 lbox
89 103 0.9500000 lbox
90 102 0.9500000 lbox
92 101 0.9500000 lbox
93 100 0.9500000 lbox
94 99 0.9500000 lbox
showpage
end
%%EOF

From `RNAfold` manpage

It also produces PostScript files with plots of the resulting secondary structure graph and a "dot plot" of the base pairing matrix. The dot plot shows a matrix of squares with area proportional to the pairing probability in the upper right half, and one square for each pair in the minimum free energy structure in the lower left half. For each pair i−j with probability p>10E−6 there is a line of the form

i j sqrt(p) ubox

in the PostScript file, so that the pair probabilities can be easily extracted.



In [7]:

    
%%bash

cat ../data/ROSE1_dp.ps | grep "^[0-9].*ubox$"









    



1 16 0.988398700 ubox
1 46 0.004186278 ubox
1 52 0.006048086 ubox
1 109 0.006282958 ubox
1 114 0.010033249 ubox
1 116 0.019380207 ubox
2 15 0.999615089 ubox
2 102 0.006665712 ubox
2 108 0.006124445 ubox
2 113 0.009015848 ubox
2 115 0.016778111 ubox
3 14 0.999864808 ubox
3 17 0.003515375 ubox
3 101 0.006674071 ubox
3 107 0.006041594 ubox
4 13 0.999522602 ubox
4 16 0.004822295 ubox
4 106 0.003943676 ubox
4 109 0.007687185 ubox
5 12 0.999914702 ubox
5 15 0.004874945 ubox
5 99 0.007129105 ubox
5 108 0.007697248 ubox
6 11 0.999524570 ubox
6 14 0.004875357 ubox
6 98 0.007140052 ubox
6 107 0.007696684 ubox
7 13 0.004057106 ubox
7 97 0.007089061 ubox
7 106 0.007502402 ubox
8 96 0.007073948 ubox
8 105 0.007181458 ubox
9 90 0.004925009 ubox
9 95 0.006456978 ubox
9 104 0.003520763 ubox
10 89 0.006106488 ubox
10 93 0.003873335 ubox
11 88 0.006365566 ubox
11 92 0.003807256 ubox
11 93 0.004710490 ubox
12 87 0.007402256 ubox
12 92 0.004917975 ubox
13 86 0.007501189 ubox
13 91 0.004789921 ubox
14 47 0.006886134 ubox
14 85 0.007575662 ubox
14 90 0.003271218 ubox
15 46 0.013083471 ubox
15 84 0.007465038 ubox
15 116 0.004516158 ubox
16 44 0.137293296 ubox
16 48 0.008186177 ubox
16 53 0.050426392 ubox
16 115 0.005074249 ubox
17 25 0.006034476 ubox
17 43 0.596518062 ubox
17 46 0.718032769 ubox
17 47 0.007504468 ubox
17 52 0.206020361 ubox
17 104 0.009791103 ubox
17 109 0.012884538 ubox
17 114 0.011564964 ubox
17 116 0.024623449 ubox
18 24 0.006147723 ubox
18 41 0.080928351 ubox
18 42 0.582695703 ubox
18 45 0.719111559 ubox
18 49 0.015238516 ubox
18 51 0.206216622 ubox
18 103 0.010282846 ubox
18 108 0.012776055 ubox
18 113 0.010923397 ubox
18 115 0.023411116 ubox
19 23 0.006118762 ubox
19 40 0.083427376 ubox
19 44 0.718551816 ubox
19 48 0.016439739 ubox
19 50 0.206085302 ubox
19 53 0.003809741 ubox
19 102 0.010366245 ubox
19 107 0.012579371 ubox
20 39 0.068192199 ubox
20 40 0.995071006 ubox
20 101 0.008563875 ubox
21 38 0.048232243 ubox
21 39 0.997505086 ubox
21 40 0.006866296 ubox
21 99 0.006479541 ubox
21 108 0.004549614 ubox
22 38 0.998666948 ubox
22 39 0.004892970 ubox
22 98 0.006479300 ubox
22 107 0.004626752 ubox
23 37 0.999827668 ubox
23 97 0.006343742 ubox
23 106 0.004616204 ubox
24 36 0.999558262 ubox
24 96 0.006235244 ubox
24 105 0.004545728 ubox
25 35 0.998130302 ubox
25 53 0.012233155 ubox
26 36 0.016594313 ubox
26 52 0.012631521 ubox
26 92 0.006268213 ubox
27 34 0.996329276 ubox
27 35 0.022878151 ubox
27 91 0.006262904 ubox
27 102 0.004660533 ubox
28 33 0.996667555 ubox
28 34 0.024661577 ubox
28 50 0.013659363 ubox
28 101 0.004661481 ubox
29 33 0.022389966 ubox
29 49 0.013697290 ubox
29 100 0.004620635 ubox
30 48 0.013693874 ubox
30 99 0.004459090 ubox
31 47 0.009611206 ubox
31 87 0.006030195 ubox
32 86 0.006335947 ubox
33 85 0.006362213 ubox
33 96 0.003484000 ubox
34 47 0.010173601 ubox
34 84 0.006288959 ubox
34 95 0.003519777 ubox
35 46 0.015684040 ubox
35 94 0.003383716 ubox
36 45 0.015751850 ubox
37 44 0.015782921 ubox
37 91 0.003450657 ubox
38 43 0.015740921 ubox
38 90 0.003438152 ubox
39 89 0.003328040 ubox
40 88 0.003313508 ubox
41 87 0.003274310 ubox
42 81 0.004244676 ubox
43 53 0.118752331 ubox
43 80 0.004398433 ubox
43 107 0.003191091 ubox
43 115 0.010248864 ubox
44 52 0.123479039 ubox
44 79 0.004162342 ubox
44 84 0.017665599 ubox
44 114 0.010726167 ubox
45 78 0.003620496 ubox
45 81 0.035318284 ubox
45 83 0.017433339 ubox
46 54 0.021466464 ubox
46 80 0.035854874 ubox
47 53 0.059737253 ubox
47 80 0.075419201 ubox
47 82 0.019325777 ubox
47 115 0.003529058 ubox
47 117 0.008554085 ubox
48 52 0.068667562 ubox
48 79 0.308544264 ubox
48 81 0.020759128 ubox
48 84 0.005296009 ubox
48 114 0.003728964 ubox
48 116 0.009762205 ubox
49 78 0.312243547 ubox
49 83 0.005232275 ubox
50 77 0.312441391 ubox
50 79 0.095257329 ubox
51 78 0.096333515 ubox
51 81 0.007031727 ubox
52 70 0.005694488 ubox
52 76 0.293488512 ubox
52 80 0.007525824 ubox
52 115 0.008230668 ubox
53 69 0.005814571 ubox
53 75 0.246143510 ubox
53 79 0.007630834 ubox
53 114 0.008210608 ubox
54 68 0.005815355 ubox
54 74 0.174182595 ubox
54 75 0.968096543 ubox
55 67 0.005813649 ubox
55 74 0.984681490 ubox
56 66 0.005820674 ubox
56 73 0.999926564 ubox
57 65 0.005817896 ubox
57 72 0.999773030 ubox
58 64 0.005809187 ubox
58 71 0.999972313 ubox
59 70 0.999960601 ubox
60 69 0.999948977 ubox
61 67 0.003614894 ubox
61 68 0.998269043 ubox
62 66 0.063179782 ubox
71 116 0.004021265 ubox
72 115 0.004023649 ubox
73 114 0.004029499 ubox
74 113 0.004021168 ubox
75 115 0.004449103 ubox
76 109 0.018779394 ubox
76 114 0.025320957 ubox
76 116 0.940534703 ubox
77 108 0.018718570 ubox
77 113 0.024492141 ubox
77 115 0.941184813 ubox
78 107 0.009382259 ubox
78 108 0.004345699 ubox
78 117 0.009580379 ubox
79 107 0.016983132 ubox
79 113 0.940095931 ubox
79 115 0.077522482 ubox
80 92 0.018170656 ubox
80 104 0.007369408 ubox
80 106 0.017068379 ubox
80 109 0.009183489 ubox
80 112 0.947874711 ubox
80 114 0.294118824 ubox
80 116 0.017130330 ubox
81 91 0.018205828 ubox
81 103 0.007860577 ubox
81 108 0.006818192 ubox
81 111 0.949089450 ubox
81 113 0.294482137 ubox
81 115 0.016409382 ubox
82 90 0.018020866 ubox
82 104 0.034137816 ubox
82 110 0.931742249 ubox
82 112 0.292939159 ubox
83 100 0.037590153 ubox
83 103 0.036750096 ubox
83 108 0.008044688 ubox
83 111 0.260049679 ubox
84 99 0.040093121 ubox
84 102 0.037193058 ubox
84 107 0.007395172 ubox
84 108 0.995544252 ubox
85 91 0.003384021 ubox
85 98 0.040108800 ubox
85 101 0.037091949 ubox
85 107 0.997109557 ubox
85 113 0.018621737 ubox
86 97 0.039357404 ubox
86 104 0.022824723 ubox
86 106 0.991749284 ubox
86 112 0.018686284 ubox
87 99 0.012909796 ubox
87 100 0.030830436 ubox
87 101 0.007139221 ubox
87 102 0.019874731 ubox
87 103 0.041762738 ubox
87 108 0.003291841 ubox
87 111 0.017858469 ubox
88 98 0.012408471 ubox
88 99 0.030956594 ubox
88 100 0.008806228 ubox
88 101 0.020630161 ubox
88 102 0.043113854 ubox
88 103 0.070450594 ubox
88 107 0.004035595 ubox
88 108 0.004256930 ubox
89 98 0.030643115 ubox
89 99 0.009066030 ubox
89 100 0.022056031 ubox
89 101 0.036302572 ubox
89 102 0.067975413 ubox
89 103 0.836896471 ubox
89 107 0.005390753 ubox
90 98 0.008707080 ubox
90 99 0.022196417 ubox
90 100 0.035806539 ubox
90 101 0.038577585 ubox
90 102 0.756679553 ubox
90 115 0.007650119 ubox
90 117 0.011910439 ubox
91 97 0.006359660 ubox
91 106 0.004476429 ubox
91 109 0.025497988 ubox
91 114 0.010369409 ubox
91 116 0.015146413 ubox
92 99 0.093227775 ubox
92 101 0.963924449 ubox
92 102 0.023998989 ubox
92 108 0.025530170 ubox
92 113 0.010368411 ubox
92 115 0.015133607 ubox
93 98 0.083480332 ubox
93 99 0.055047198 ubox
93 100 0.969560715 ubox
93 101 0.018014740 ubox
93 102 0.007671418 ubox
93 103 0.182213229 ubox
93 107 0.015888711 ubox
94 98 0.061578711 ubox
94 99 0.965880698 ubox
94 101 0.012504460 ubox
94 102 0.193292539 ubox
95 101 0.193501122 ubox
95 102 0.004554297 ubox
96 100 0.174327583 ubox
96 101 0.004449040 ubox
98 110 0.003230757 ubox
99 109 0.003226014 ubox
100 110 0.003269619 ubox
101 109 0.003822609 ubox
101 110 0.005433951 ubox
102 109 0.005653989 ubox
106 115 0.007642579 ubox
106 117 0.017802702 ubox
107 114 0.011088268 ubox
107 116 0.022419345 ubox
108 114 0.031753481 ubox
108 116 0.044261073 ubox
109 113 0.031642938 ubox
109 115 0.044659139 ubox



In [16]:

    
%%bash

awk '/^>/' ../data/rose.fa | head -1









    



>ROSE1



In [ ]:

    
%%bash


# Runs RNAfold for the given RNA sequence over the range of temperatures,
# extracts base pairing probabilities and saves them in .txt files 
# for later analysis

# PostScript file generated by RNAfold ends with this
psext="_dp.ps"

# FASTA file with RNA sequence, rna.fa by default
rna_fa=${1:-rna.fa}

# Temperature interval limits
T1=${2:-37}
T2=${3:-43}

# Check the input file exists
if [[ ! -f $rna_fa ]]
then
    echo "Could not find $rna_fa ... Exiting."
    exit 1
fi

# Get the base_name either from the fasta file or the filename
base_name=`awk '/^>/' $rna_fa | head -1`
if [[ -z "$base_name" ]]
then
    base_name="${rna_fa%.*}"
else
    base_name=${base_name##>}
fi
                
# Iterate over the T range and save probabilities to .txt file
for T in $(seq $T1 $T2)
do
    echo "Running RNAfold for Temp=$T ..."
    RNAfold -p -d2 --noPS --noLP -T $T < $rna_fa
    tmpf=`ls | grep _dp.ps`
    grep "^[0-9].*ubox$" $tmpf > ${base_name}_${T}.txt
done
               
# Cleanup
rm ${base_name}_dp.ps



In [ ]:

From RNAfold manpage

From `RNAfold` manpage