Mapping Reads to a Reference Genome

Shell Variables

[1]:
# Source the config script
source bioinf_intro_config.sh

ls $CUROUT
count_out  igv     qc_output  stuff_for_igv_shorter_intron.tgz  trimmed_fastqs
genome     myinfo  star_out   stuff_for_igv.tgz

Mapping with STAR

[2]:
STAR \
    --runMode alignReads \
    --twopassMode None \
    --genomeDir $GENOME_DIR \
    --readFilesIn $TRIMMED/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz \
    --readFilesCommand gunzip -c \
    --outFileNamePrefix ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_ \
    --quantMode GeneCounts \
    --outSAMtype BAM Unsorted \
    --outSAMunmapped Within \
    --runThreadN 2
Jun 26 15:10:01 ..... started STAR run
Jun 26 15:10:01 ..... loading genome
Jun 26 15:10:05 ..... started mapping
Jun 26 15:11:40 ..... finished successfully

We will start with these parameters, but there is an extensive list of command line options detailed in the STAR Manual, it is a good idea to read through and try to understand all of them. We will discuss some more later.

  • –runMode alignReads : map reads

  • –twopassMode : Run one pass or two? If two-pass mode is on, STAR tries to discover novel junctions, then reruns mapping with these added to the annotation

  • –genomeDir : directory containing the genome index

  • –readFilesIn : input FASTQ

  • –readFilesCommand gunzip -c : use “gunzip -c” to uncompress FASTQ on-the-fly, since it is gzipped

  • –outFileNamePrefix : prefix (and path) to use for all output files

  • –quantMode GeneCounts : output a table of read counts per gene

  • –outSAMtype BAM Unsorted : output an unsort BAM file

  • –outSAMunmapped Within : included unmapped reads in the BAM file

  • –runThreadN : tells STAR to run using multiple cores. I am using it so we don’t have to wait too long for this to run during class. It is OK to use multiple cores, but before you do this you should be sure that the server is not busy, and even then you should use a reasonable number of cores. Abusing multi-threading is inconsiderate of other users and could crash the server.

STAR Output

So what happened? Let’s take a look …

[3]:
ls ${STAR_OUT}
21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L001_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.out
21_2019_P_M1_S21_L001_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L001_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L001_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_Aligned.out.bam
21_2019_P_M1_S21_L002_R1_Log.final.out
21_2019_P_M1_S21_L002_R1_Log.out
21_2019_P_M1_S21_L002_R1_Log.progress.out
21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam
21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai
21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.out
21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out
21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab
21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab
21_2019_P_M1_S21_L002_R1_SJ.out.tab
genome_Log.out
[4]:
head ${STAR_OUT}/21_2019_P_M1_S21_L002_R1*
==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Aligned.out.bam <==
�BC���Mo�@���C
$�@�H�X!���_�S�T�JM��M���^����

�?�g� =��r�#�q�R� �#�M�AX�qLlӪVq,�1gy��^�4�Z{p_(���Z�2b2V2O4�� �)�

�Lk� ���d&amp;0�.��,`��0������v��u�5�p��#�ڬb qm�e-��Ř=�

6:�Z�W���1?��a5��]�����z� /

q=� �b�����D {:0�”��r���D �q�t?�^�&lt;1�a����

�2ө������ �f����Y�eF��A*wS��b�����NkS���F�ۀ���@5��g7��@��q&amp;�i

B��=�d&gt;V/�x?��ᳳN�� �U���TƠ���)ח:�BD’���

�gB�ިO�.��i���ł���{Z��D

��

�5YY�y��ϗ��^����U]�j&amp;��s�9������40=b�;��

#�ĮMEKGC QT;a.pF1�AŏN����cl�J�a���Y�0���Uu�{�ǩ�Y+����ԩ�{������~��g��~����0{��^ ������Twp:{�Y�͏�����8|F�����s^b��x�|�����x��������G�fy^�Ɩ���&amp;��8I��[$u�DMEm��ud�q՝������’��������Vw�T=�Λ����W�����c����`z�n�ӒC�� �(lӠ��1���aS�eX��d�i�*�L?�S#&lt;��i̙�j���C����)tZ��{�t�0�i�,HB�%a��IPU�f6.m��”�y$n0P�H��P:���Q�^���ʕk�r%��]�~���kׯ�l�ӣ��_x����I�櫽� &amp;DEaY�Q��^�E �MY��j��$�q��咇��lwq�N�

�[���e���O�c��.&gt;�0�VLJ

�;4��06a֤Un�, �866��&amp;/�0,�”/��o��&amp;�#�bl��U����܍ۿ�k����k :e#D��l�� ?�&lt;D�h

�V��� �&amp;,�6o”�D��G�Viۤ��”�aRY/ ����

� ��{���8ݣ#HN@�p ���W`�&amp;��&lt;&gt;�@�� 2*

��J���4 3 ��&amp;(�m�ب��U)׃X@�׮�!�h�t�N����

8tW�z1��̛�9���W&lt;t+��}

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.final.out &lt;==

Started job on | Jun 26 15:10:01

Started mapping on | Jun 26 15:10:05

Finished on | Jun 26 15:11:40

Mapping speed, Million of reads per hour | 92.30

Number of input reads | 2435761

Average input read length | 75

UNIQUE READS:

Uniquely mapped reads number | 2354695

Uniquely mapped reads % | 96.67%

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.out &lt;== STAR version=STAR_2.5.2b STAR compilation time,server,dir=&lt;not set in Debian&gt; ##### DEFAULT parameters: versionSTAR 20201 versionGenome 20101 20200 parametersFiles - sysShell - runMode alignReads runThreadN 1 runDirPerm User_RWX

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.progress.out &lt;==
Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped

M/hr number length unique length MMrate multi multi+ MM short other

Jun 26 15:11:06 94.5 1601279 75 96.7% 75.4 0.2% 1.4% 0.1% 0.0% 1.8% 0.0% ALL DONE!

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab &lt;== N_unmapped 46060 46060 46060 N_multimapping 35006 35006 35006 N_noFeature 12466 2145291 18783 N_ambiguous 203327 820 316 CNAG_04548 0 0 0 CNAG_07303 0 0 0 CNAG_07304 6 0 6 CNAG_00001 0 0 0 CNAG_07305 1 0 1 CNAG_00002 51 0 51

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam &lt;== �F��T����i����Di���(��{y�ׇ��L�v33���S�[��f�nJ ^B�;�ݝ�ɬ}A�я�I��z�#������� Ba����̓�Fp��u��,�:��F�V�i�KD,ye!|f�j�Xm�)e�f�`�h3R”:%ĥ����Հ9�ClR0�m�&gt;��V8�u�ٔ�n�`

�c.��|���#��c�&lt;� ;��v �ޡK;�pϳ�KKG� ׳�}pfS��* �ece� �7:F���

`�E:�KIfZC%���XX�i��ʻY,�Ck+Iя!���z25��d��q�j�%�1�E�8�;�EID1aQ; H��*�4cR�z�`�LxN��L��X���D

$�ӱn�I’rji�c��x?ٻ��rg��&gt;_�0ͤv�&lt;*V �����b#P(T�[�o�I?���i?I���A���w�: ���`h�nn��O}��$�x� �L’B�E��M��&lt;��V*ԯ�3&amp;9����g%��M�˙a9=��咣r~JTF.3 �ҩ̜�������Q�=��Wn!��������kJm�R���}P5�kF+��8 �A?�΁rћowAy��V&gt;���C�n���׏f�$���{ZR�����54��H��N�&amp;� �Ъ�.�BCe=�}{�,[VV�{�0�h���:�z�9Ֆڀ@�_�s�ny���+�0*��#� ����11��Ҵ&lt;�����

��u�1b4����׮]���Uu�ĩ{�U��O�W��k����{���

���F��{�����{Q��������b?�W��1u�{5}����X.��b�����~5��K�QEdyz�� &amp;�)+�I0 �sM��y���o&lt;�_ �E�&gt;�ǟ���_y)�+/&amp;_��{���2��h�݀��#�}]!-���߮���b�} ����ir�ϓ�)���i��t��i��8 G�N� ��{�@z; } }����,K�!-i{�W��j7�4ˣl�f�8�&amp;�I�f��li:�r�F�l�N�Cz��nKΓXY�o�~R=��cI }²�.@R��i�V��q��n������HNI���p�M�I�g�)���&lt;b�`��ݩy�bzIi��ӟ&amp;&amp;+y_�G��� mi�ޭ�� z��O�,��s�y��Dc bBz����������8@��* &gt;$I���’-3OP&amp;x

����)�7 ���[���?��6����0?CA�I2���l�AI��4:EA ‘Y;S�Z�-�B�����c�!j�’�@-��������j��DM�u

�� OF)~ڢj��S���a�ܬ�#� ��Z4�mG�qL�i���( �b�J_�c ���MP���-��������q�8,�����Ί�*o��`����i���������v =S��(, �{�m�w��&gt;�7”�=}��%��5�N��f�[�|�nb̻ ̻S8���

#�����ڕ���a=z���m�~t(XF�T�a=7�����0�a.����%���f��%2؎=u㠒x��О`-�[���0L=P�j��u� w������Z-�����򸜯��5n&gt;m���&amp;�����G”Ӧ^�ug+��

�|�_ÉZ�a”��9}�8�I’ERde���FPx콒4+2|�K�-�&gt;=��� hW�NԿ�*��o=��n���EO �W��T�

�6�) �����u�&lt;[����c �e� ` �������}ᓠ���?^A��ӂ·��#z`c�O�K�`.����F�d*U�U�Gk��P

#x&gt;��D���OX�GL�:��:�d � �����ކ=PHo&gt;  ��Wsw �ԏ��������~+

��đև�R�}������R�29�C�{�g�G���Ef��/������

5p�TQ��kS=��55�b�[�ׅ� ��2�z�II ?l*

�{��!����Zn;t‹�A�DC1�S�e

��s��� x��”=&lt;|�K�z�|C�� P��#�{�F}��ٜ�#&lt;L��Q�g�� ��v�v{6���O�a:� ��br�� C*��9��������7yz�g@���s���at]wt� ��� �� �fE���&gt;UHlS��+�[@z��$��ǭ���S��_��CO@�i��I���� �F�����,I�I��)&gt;’Y N�&lt;���$/�@3J��S�.� P9ƾxca�C&amp;`џb/������߻˨�����8(��r�H�|]�~I���8,ð ���9��E8��0��cj��z�c4��l�#���Ǐ�t��-X&lt;/�T� � ���Z_����|u6�`Ԓ���~s0ͪ����#O�� �r[�?$�t%K����WaB��G�:�n����+�d� �3�.CT��za�”����C�������7*�����1��_/��J����[[��(���7���Y�~�����z6�`�#J1���G�{VLZ�N߭�4&gt;F��3Ջ����|�; �qn�v

B#�x� �:��A|յ����

���[,w{c�u��� y*�� 4�T�U��I��

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai &lt;== �JM����(�c�1+���B,�Έ+-rpT/+�pT/t�a0�rU�0`�0ϱ

���^��N��7^&gt; n�Έ+-rpT/o+�pT/t�a0l�c�1+���B,rsp�1�Y�;5z���N&lt;qMj=~^�LCZ �`D��j�EO[�F}f �]&gt;�cLCN4yU�0sp�1��M6��av9O7���;���N&lt;qMj=f �]&gt;X,LC^�LCZ �`D�j�Ex��av9[���;s�Y�;5��M6�nE_pG�ڣ�I�!�jvM��JQ���JQ�.K�TP�2[�FnE_pG�ڣ�I�(L�L��|KV�~:�W��:�WVq��X�Vq��X���Z����Zk.h�^�&gt;�T�k�Y��l�W�

`+xcQ�*L�L!�jvM�.K�T�|KV�~:�W�:�WRI���_W� `I’�e��բf�k.h�^I���_���բf��g�+xcI’�eJ���Z���Z��g� ��g&lt;� o�� o󉱎x����x50�d�05�d�����h�

i�^� i��T�k�&lt;&amp;xPqJʹrS� ��g���h

�5��{yl~5T�̂U{�nz�m6Q{��2`|馷}Zpl~{yl~5T�̂50�d���m6Q{��2`|���5�,lYÎ�׉��@ҋ��q�t��B����B��05�d��0&gt;���0&gt;�� ��?���J΢�JS�j��R�f”N��N��7^&gt;

���� ����� 1��!� x�&amp;��J�q:1���ޚ������ƹ�d�d���=x���`�ի��:1�G�g;���g;��=ƹ���,�(�JAD��X��xٜ�J�q:1��:1�G�g;��g;��6ƹ��ƹ�d�d��^e�$�Y^e�$�=x���`�ի8�ի��&amp;T��,�(�JAD���d�Y�Z�d�Y��n���������ζs��8�ի��&amp;T���n��������J���ζs�!���� @�%��r�h�����������y�/�w����#W��k� :+�j ��դ �E�

�6’�yz’��8�’���q(����(�@��)2��*4b�1+���B,�Έ+-rpT/t�a0�rU�0�p�1�Y�;5��M6���6���7XE��8��av9[���;���N&lt;qMj=����=��]&gt;X,LCZ

�`DFf�EO[�FnE_pG�ڣ�I��9’J�ğjJ�*�Js�`K��w�K�(L�L!�jvM��JQ�.K�TL��U�|KV�~:�WVq��X���Zk.h�^I���_�

`+xcI’�evR�+f��բf��g���h�5�,lYÎ�5������@ҋ�`�t��B��o���gpZ%xPqJʹr���t�S�v��*wy�7x󉱎x{�nz�m6Q{��2`|ɘ�}Zpl~5T�̂�&lt;g��50�d��0&gt;��

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out &lt;==

Started job on | Jun 24 09:50:39

Started mapping on | Jun 24 09:50:39

Finished on | Jun 24 09:51:17

Mapping speed, Million of reads per hour | 230.76

Number of input reads | 2435761

Average input read length | 75

UNIQUE READS:

Uniquely mapped reads number | 2355964

Uniquely mapped reads % | 96.72%

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.out &lt;== STAR version=STAR_2.5.2b STAR compilation time,server,dir=&lt;not set in Debian&gt; ##### DEFAULT parameters: versionSTAR 20201 versionGenome 20101 20200 parametersFiles - sysShell - runMode alignReads runThreadN 1 runDirPerm User_RWX

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out &lt;==
Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped

M/hr number length unique length MMrate multi multi+ MM short other

ALL DONE!

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab &lt;== N_unmapped 46059 46059 46059 N_multimapping 33738 33738 33738 N_noFeature 12438 2146527 18736 N_ambiguous 203237 786 143 CNAG_04548 0 0 0 CNAG_07303 0 0 0 CNAG_07304 6 0 6 CNAG_00001 0 0 0 CNAG_07305 1 0 1 CNAG_00002 51 0 51

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab &lt;== 1 9247 9310 2 2 1 6 4 35 1 10455 10506 1 3 0 7 8 24 1 10609 10694 1 1 0 0 2 21 1 11574 11622 2 2 0 0 1 21 1 11574 11637 2 2 0 0 2 32 1 12904 12991 1 1 1 4 0 37 1 13077 13136 1 1 1 10 0 37 1 13256 13310 1 1 1 7 0 37 1 13480 13530 1 1 1 5 0 33 1 13480 13534 1 1 0 1 0 28

==&gt; /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_SJ.out.tab &lt;== 1 9247 9310 2 2 1 6 4 35 1 10455 10506 1 3 0 7 7 24 1 10609 10694 1 1 0 0 2 21 1 11574 11622 2 2 0 0 1 21 1 11574 11637 2 2 0 0 2 32 1 12904 12991 1 1 1 4 0 37 1 13077 13136 1 1 1 10 0 37 1 13256 13310 1 1 1 7 0 37 1 13480 13530 1 1 1 5 0 33 1 13480 13534 1 1 0 1 0 28 </pre>

$�@�H�Xtextbackslash{}!���_�S�T�JM��M���^{}����

�?�g� =��r�#�q�R� �#�M�AX�qLlӪVq,�1gy��^{}�4�Z{p_(���Z�2b2V2O4�� �)�

�Lk� ���d&0�.��,`��0������v��u�5�p��#�ڬb qm�e-��Ř=�

6:�Z�W���1?��a5��]�����z� /

q=� �b�����D {:0�”��r���D �q�t?�^{}�<1�a����

�2ө������ �f����Y�eF��A*wS��b�����NkS���F�ۀ���@5��g7��@��q&�i

B��=�d>V/�x?��ᳳN�� �U���TƠ���)ח:�BD’���

�gB�ިO�.��i���ł���{Z��D

��

�5YY�y��ϗ��^{}����U]�j&��s�9������40=b�;��

#�ĮMEKGC QT;atextbackslash{}.pF1�AŏN����cl�J�a���Y�0���Uu�{�ǩ�Y+����ԩ�{������textasciitilde{}��g��textasciitilde{}����0{��^{} ������Twp:{�Y�͏�����8|F�����s^{}b��x�|�����x��������G�fy^{}�Ɩ���&��8I��[$u�DMEm��ud�q՝������’��������Vw�T=�Λ����W�����c����`z�n�ӒC�� �(lӠ��1���aS�eX��d�i�*�L?�S#<��i̙�j���C����)tZ��{�t�0�i�,HB�%a��IPU�f6.m��”�y$n0P�H��P:���Q�^{}���ʕk�r%��]�textasciitilde{}���kׯ�l�ӣ��_x����I�櫽� &DEaY�Q��^{}�E �MY��j��$�q��咇��lwq�N�

�[���e���O�c��.>�0�VLJ

�;4��06a֤Un�, �866��&/�0,�”/��o��&�#�bl��U����܍ۿ�k����k :e#D��l�� ?�<D�h

�V��� �&,�6o”�D��G�Viۤ��”�aRY/ ����

� ��{���8ݣ#HN@�p ���W`�&��<>�@�� 2*

��J���4 3 ��&(�m�ب��U)׃X@�׮�!�h�t�N����

8tW�z1��̛�9���W<t+��}

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.final.out <==

Started job on | Jun 26 15:10:01

Started mapping on | Jun 26 15:10:05

Finished on | Jun 26 15:11:40

Mapping speed, Million of reads per hour | 92.30

Number of input reads | 2435761

Average input read length | 75

UNIQUE READS:

Uniquely mapped reads number | 2354695

Uniquely mapped reads % | 96.67%

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.out <== STAR version=STAR_2.5.2b STAR compilation time,server,dir=<not set in Debian> ##### DEFAULT parameters: versionSTAR 20201 versionGenome 20101 20200 parametersFiles - sysShell - runMode alignReads runThreadN 1 runDirPerm User_RWX

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.progress.out <==
Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped

M/hr number length unique length MMrate multi multi+ MM short other

Jun 26 15:11:06 94.5 1601279 75 96.7% 75.4 0.2% 1.4% 0.1% 0.0% 1.8% 0.0% ALL DONE!

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab <== N_unmapped 46060 46060 46060 N_multimapping 35006 35006 35006 N_noFeature 12466 2145291 18783 N_ambiguous 203327 820 316 CNAG_04548 0 0 0 CNAG_07303 0 0 0 CNAG_07304 6 0 6 CNAG_00001 0 0 0 CNAG_07305 1 0 1 CNAG_00002 51 0 51

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam <== �F��T��textbackslash{}��i����Di���(��{y�ׇ��L�v33���S�[��f�nJ ^{}B�;�ݝ�ɬ}A�я�I��z�#������� Ba����̓�Fp��u��,�:��F�V�i�KD,ye!|f�j�Xm�)e�f�`�h3R”:%ĥ����Հ9�ClR0�m�>��V8�u�ٔ�n�`

�c.��|���#��c�<� ;��v �ޡK;�pϳ�KKG� ׳�}pfS��* �ece� �7:F���

`�E:�KIfZC%���XX�i��ʻY,�Ck+Iя!���z25��d��q�j�%�1�E�8�;�EID1aQ; H��*�4cR�z�`�LxN��L��X���D

$�ӱn�I’rji�c��x?ٻ��rg��>_�0ͤv�<*V �����b#P(T�[�o�I?���i?I���A���w�: ���`h�nn��O}��$�x� �L’B�E��M��<��V*ԯ�textbackslash{}3&9����g%��M�˙a9=��咣rtextasciitilde{}JTF.3 �textbackslash{}ҩ̜�������Q�=��Wn!��������kJm�R���}P5�kF+��8 �A?�΁rћowAy��V>���C�n���׏f�$���{ZR�����54��H��N�&� �Ъ�.�BCe=�}{�,[VV�{�0�h���:�z�9Ֆڀ@�_�s�ny���+�0*��#� ����11��Ҵ<�����

��u�1b4����׮]���Uu�ĩ{�U��O�W��k����{���

���F��{�����{Q��������b?�W��1u�{5}����X.��b�����textasciitilde{}5��K�QEdyz�� &�)+�I0 �sM��y���o<�_ �E�>�ǟ���_y)�+/&_��{���2��h�݀��#�}]!-���߮���b�} ����ir�ϓ�)���i��t��i��8 G�N� ��{�@z; } }����,K�!-i{textbackslash{}�W��j7�4ˣl�f�8�&�I�f��li:�r�F�l�N�Cz��nKΓXY�o�textasciitilde{}R=��cI }²�.@R��i�V��q��n������HNI���p�M�I�g�)���<b�`��ݩy�bzIi��ӟ&&+y_�G��� mi�ޭ�� z��O�,��s�y��Dc bBz����������8@��* >$I���’-3OP&x

����)�7 ���[���?��6����0?CA�I2���l�AI��4:EA ‘Y;S�Z�-�B�����c�!j�’�@-��������j��DM�u

�� OF)textasciitilde{}ڢj��S���a�ܬ�#� ��Z4�mG�qL�i���( �b�J_�c ���MP���-��������q�8,�����Ί�*o��`����i���������v =S��(, �{�m�w��>�7”�=}��%��5�N��f�[�|�nb̻ ̻S8���

#�����ڕ���a=z���m�textasciitilde{}t(XF�T�a=7�����0�a.����%���f��%2؎=u㠒x��О`-�[���0L=P�j��u� w������Z-�����򸜯��5n>m���&�����G”Ӧ^{}�ug+��

�|�_ÉZ�a”��9}�8�I’ERde���FPx콒4+2|�K�-�>=��� hW�NԿ�*��o=��n���EO �W��T�

�6�) �����u�<[����c �e� ` �������}ᓠ���?^{}A��ӂ·textbackslash{}��#z`c�O�K�`.����textbackslash{}F�d*textbackslash{}U�U�Gk��P

#x>��D���OX�GL�:��:�d � �����ކ=PHo>  ��Wsw �ԏ�����textbackslash{} ���textasciitilde{}+

��đև�R�}�����textbackslash{}�R�29�C�{�g�G���Ef��/������

5p�TQ��kS=��55�b�[�ׅ� ��2�z�II ?l*

�{��!����Zn;t‹�A�DC1�S�e

��s��� x��”=<|�K�z�|C�� P��\#�\{�F\}��ٜ�\#<L��Q�g�� ��v�v\{6���O�a:� ��br�� C*��9��������7yz�g@���s���at]wt� ��� �� �fE���>UHlS��+�[@z��\$��ǭ���S��\_��CO@�i��I\textbackslash{}���� �F�����,I�I��)>'Y N�<���\$/�@3J��S�.� P9ƾxca�C\&`џb/������߻\˨�����8(��r�H�|]�textasciitilde{}I���8,ð ���9��E8��0��cj��z�c4��l�#���Ǐ�t��-X</�T� � ���Z_����|u6�`Ԓ���textasciitilde{}s0ͪ����#O�� �r[�?$�t%K����WaB��G�:�n����+�d� �3�.CT��za�”����C�������7*�����1��_/��J����[[��(���7���Y�textasciitilde{}�����z6�`�#J1���G�{VLZ�N߭�4>F��3Ջ����|�; �qn�v

B#�x� �:��A|յ����

���[,w{c�u��� y*�� 4�T�U��I��

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai <== �JM����(�c�1+���B,�Έ+-rpT/+�pT/t�a0�rU�0`�0ϱ

���^{}��N��7^{}> n�Έ+-rpT/o+�pT/t�a0l�c�1+���B,rsp�1�Y�;5z���N<qMj=textasciitilde{}^{}�LCZ �`D��j�EO[�F}f �]>�cLCN4yU�0sp�1��M6��av9O7���;���N<qMj=f �]>X,LC^{}�LCZ �`D�j�Ex��av9[���;s�Y�;5��M6�nE_pG�ڣ�I�!�jvM��JQ���JQ�.K�TP�2[�FnE_pG�ڣ�I�(L�L��|KV�textasciitilde{}:�W��:�WVq��X�Vq��X���Z����Zk.h�^{}�>�T�k�Y��l�W�

textbackslash{}`+xcQ�*L�L!�jvM�.K�T�|KV�textasciitilde{}:�W�:�WRI���_W� textbackslash{}`I’�e��բf�k.h�^{}I���_���բf��g�+xcI’�eJ���Z���Z��g� ��g<� o�� o󉱎x����x50�d�05�d�����h�

i�^{}� i��T�k�<&xPqJʹtextbackslash{}rS� ��g���h

�5��{yltextasciitilde{}5T�̂U{�nz�m6Q{��2`|馷}Zpltextasciitilde{}{yltextasciitilde{}5T�̂50�d���m6Q{��2`|���5�,lYÎ�׉��@ҋ��q�t��B����B��05�d��0>���0>�� ��?���J΢�JS�j��R�f”N��N��7^{}>

���� ����� 1��!� x�&��J�q:1���ޚ������ƹ�d�d���=x���`�ի��:1�G�g;���g;��=ƹ���,�(�JAD��X��xٜ�J�q:1��:1�G�g;��g;��6ƹ��ƹ�d�d��^{}e�$�Y^{}e�$�=x���`�ի8�ի��&T��,�(�JAD���d�Y�Z�d�Y��n���������ζs��8�ի��&T���n��������J���ζs�!���� @�%��r�h�����������y�/�w����#W��k� :+�j ��դ �E�

�6’�yz’��8�’���q(����(�@��)2��*4b�1+���B,�Έ+-rpT/t�a0�rU�0�p�1�Y�;5��M6���6���7XE��8��av9[���;���N<qMj=����=��]>X,LCZ

�`DFf�EO[�FnE_pG�ڣ�I��9’J�ğjJ�*�Js�`K��w�K�(L�L!�jvM��JQ�.K�TL��U�|KV�textasciitilde{}:�WVq��X���Zk.h�^{}I���_�

textbackslash{}`+xcI’�evR�+f��բf��g���h�5�,lYÎ�5������@ҋ�`�t��B��o���gpZ%xPqJʹtextbackslash{}r���t�S�v��*wy�7x󉱎x{�nz�m6Q{��2`|ɘ�}Zpltextasciitilde{}5T�̂�<g��50�d��0>��

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out <==

Started job on | Jun 24 09:50:39

Started mapping on | Jun 24 09:50:39

Finished on | Jun 24 09:51:17

Mapping speed, Million of reads per hour | 230.76

Number of input reads | 2435761

Average input read length | 75

UNIQUE READS:

Uniquely mapped reads number | 2355964

Uniquely mapped reads % | 96.72%

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.out <== STAR version=STAR_2.5.2b STAR compilation time,server,dir=<not set in Debian> ##### DEFAULT parameters: versionSTAR 20201 versionGenome 20101 20200 parametersFiles - sysShell - runMode alignReads runThreadN 1 runDirPerm User_RWX

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out <==
Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped

M/hr number length unique length MMrate multi multi+ MM short other

ALL DONE!

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab <== N_unmapped 46059 46059 46059 N_multimapping 33738 33738 33738 N_noFeature 12438 2146527 18736 N_ambiguous 203237 786 143 CNAG_04548 0 0 0 CNAG_07303 0 0 0 CNAG_07304 6 0 6 CNAG_00001 0 0 0 CNAG_07305 1 0 1 CNAG_00002 51 0 51

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab <== 1 9247 9310 2 2 1 6 4 35 1 10455 10506 1 3 0 7 8 24 1 10609 10694 1 1 0 0 2 21 1 11574 11622 2 2 0 0 1 21 1 11574 11637 2 2 0 0 2 32 1 12904 12991 1 1 1 4 0 37 1 13077 13136 1 1 1 10 0 37 1 13256 13310 1 1 1 7 0 37 1 13480 13530 1 1 1 5 0 33 1 13480 13534 1 1 0 1 0 28

==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_SJ.out.tab <== 1 9247 9310 2 2 1 6 4 35 1 10455 10506 1 3 0 7 7 24 1 10609 10694 1 1 0 0 2 21 1 11574 11622 2 2 0 0 1 21 1 11574 11637 2 2 0 0 2 32 1 12904 12991 1 1 1 4 0 37 1 13077 13136 1 1 1 10 0 37 1 13256 13310 1 1 1 7 0 37 1 13480 13530 1 1 1 5 0 33 1 13480 13534 1 1 0 1 0 28 end{sphinxVerbatim} } % The following relax is needed to avoid problems with adjacent ANSI % cells and some other stuff (e.g. bullet lists) following ANSI cells. % See https://github.com/sphinx-doc/sphinx/issues/3594 relax

STAR generates several files for each FASTQ: * Log.out : lots of details of the run, including all parameters used * Log.final.out : Important summary statistics * ReadsPerGene.out.tab : Count table, the main thing we are interested in * SJ.out.tab : All splice junctions, including ones from the GTF and novel junctions discovered by STAR * Log.progress.out: run statistics updated during run, not so interesting at the end

Let’s take a closer look at Log.final.out

[5]:
cat ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_Log.final.out
                                 Started job on |       Jun 26 15:10:01
                             Started mapping on |       Jun 26 15:10:05
                                    Finished on |       Jun 26 15:11:40
       Mapping speed, Million of reads per hour |       92.30

                          Number of input reads |       2435761
                      Average input read length |       75
                                    UNIQUE READS:
                   Uniquely mapped reads number |       2354695
                        Uniquely mapped reads % |       96.67%
                          Average mapped length |       75.36
                       Number of splices: Total |       610596
            Number of splices: Annotated (sjdb) |       596312
                       Number of splices: GT/AG |       603894
                       Number of splices: GC/AG |       5846
                       Number of splices: AT/AC |       112
               Number of splices: Non-canonical |       744
                      Mismatch rate per base, % |       0.16%
                         Deletion rate per base |       0.01%
                        Deletion average length |       1.35
                        Insertion rate per base |       0.00%
                       Insertion average length |       1.11
                             MULTI-MAPPING READS:
        Number of reads mapped to multiple loci |       35006
             % of reads mapped to multiple loci |       1.44%
        Number of reads mapped to too many loci |       1263
             % of reads mapped to too many loci |       0.05%
                                  UNMAPPED READS:
       % of reads unmapped: too many mismatches |       0.00%
                 % of reads unmapped: too short |       1.83%
                     % of reads unmapped: other |       0.00%
                                  CHIMERIC READS:
                       Number of chimeric reads |       0
                            % of chimeric reads |       0.00%

Sanity Check: Number of Reads

Are the number of input reads what we expect? Let’s look at how many reads are in the input FASTQ

[6]:
zcat $TRIMMED/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz | awk '{s++}END{print s/4}'
2435761

Sanity Check: Unmapped Reads

It is always a good idea to examine a sample of unmapped reads to figure out what they are. The easiest way to do this is with BLAST. In the past I have discovered that an experiment was contaminated with a different species by BLASTing unmapped reads. In that case there were a large number of unmapped reads, which raised my suspicions.

Even with a high rate of mapped reads, it is worth spending a few minutes to check them out. The simplest thing to do is use samtools to generate a FASTA from the unmapped.bam, grab a few of these sequences, and then BLAST them against the nr database

[7]:
samtools fasta -f 4 ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_Aligned.out.bam | head -n20
>NB501800:327:HF27FBGXB:2:11202:18025:8102
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:20718:8137
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTAAAAGTGGGGGTTGTTTTTTATTTTTTTGTAGATTTAAAAA
>NB501800:327:HF27FBGXB:2:11202:9226:8227
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:17164:8437
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:25989:8488
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTAACTCGTATGCCGTCTTATGCTTGAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:11951:8497
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:11505:8555
ATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:21516:8563
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAATAAG
>NB501800:327:HF27FBGXB:2:11202:8434:8590
GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA
>NB501800:327:HF27FBGXB:2:11202:12081:8652
GATCGAAAGACCACACGTCTGAACACCAGTCACAGCTCAATATCACGTATGAAGTATTATGCATTAAAAAAAAA
[bam2fq_mainloop] Error writing to FASTx files.: Broken pipe
[M::bam2fq_mainloop] discarded 0 singletons
[M::bam2fq_mainloop] processed 71 reads
samtools bam2fq: Error closing STDOUT: Broken pipe

MultiQC

MultiQC also works with STAR reports, so let’s try it!

[8]:
multiqc ${STAR_OUT} --outdir ${STAR_OUT}
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:45: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
  configs = yaml.load(f)
/home/jovyan/.local/lib/python3.5/site-packages/multiqc/utils/config.py:51: YAMLLoadWarning: calling yaml.load() without Loader=... is deprecated, as the default Loader is unsafe. Please read https://msg.pyyaml.org/load for full details.
  sp = yaml.load(f)
[INFO   ]         multiqc : This is MultiQC v1.7
[INFO   ]         multiqc : Template    : default
[INFO   ]         multiqc : Searching '/home/jovyan/work/scratch/bioinf_intro/star_out'
Searching 21 files..  [####################################]  100%
[INFO   ]            star : Found 3 reports and 3 gene count files
[INFO   ]         multiqc : Compressing plot data
[INFO   ]         multiqc : Report      : ../../scratch/bioinf_intro/star_out/multiqc_report.html
[INFO   ]         multiqc : Data        : ../../scratch/bioinf_intro/star_out/multiqc_data
[INFO   ]         multiqc : MultiQC complete

Once multiqc is done running we can view the results by finding the output in the Jupyter browser, it should be in a file named multiqc_report.html in :

[9]:
echo ${STAR_OUT}
/home/jovyan/work/scratch/bioinf_intro/star_out