{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Mapping Reads to a Reference Genome\n", "\n", "\n", "## Shell Variables" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\u001b[0m\u001b[01;34mcount_out\u001b[0m \u001b[01;34migv\u001b[0m \u001b[01;34mqc_output\u001b[0m \u001b[01;31mstuff_for_igv_shorter_intron.tgz\u001b[0m \u001b[01;34mtrimmed_fastqs\u001b[0m\n", "\u001b[01;34mgenome\u001b[0m \u001b[01;34mmyinfo\u001b[0m \u001b[01;34mstar_out\u001b[0m \u001b[01;31mstuff_for_igv.tgz\u001b[0m\n" ] } ], "source": [ "# Source the config script\n", "source bioinf_intro_config.sh\n", "\n", "ls $CUROUT" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mapping with STAR\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Jun 26 15:10:01 ..... started STAR run\n", "Jun 26 15:10:01 ..... loading genome\n", "Jun 26 15:10:05 ..... started mapping\n", "Jun 26 15:11:40 ..... finished successfully\n" ] } ], "source": [ "STAR \\\n", " --runMode alignReads \\\n", " --twopassMode None \\\n", " --genomeDir $GENOME_DIR \\\n", " --readFilesIn $TRIMMED/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz \\\n", " --readFilesCommand gunzip -c \\\n", " --outFileNamePrefix ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_ \\\n", " --quantMode GeneCounts \\\n", " --outSAMtype BAM Unsorted \\\n", " --outSAMunmapped Within \\\n", " --runThreadN 2" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will start with these parameters, but there is an extensive list of command line options detailed in the [STAR Manual](https://github.com/alexdobin/STAR/blob/master/doc/STARmanual.pdf), it is a good idea to read through and try to understand all of them. We will discuss some more later.\n", "\n", "* --runMode alignReads : map reads \n", "* --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\n", "* --genomeDir : directory containing the genome index\n", "* --readFilesIn : input FASTQ\n", "* --readFilesCommand gunzip -c : use \"gunzip -c\" to uncompress FASTQ on-the-fly, since it is gzipped\n", "* --outFileNamePrefix : prefix (and path) to use for all output files\n", "* --quantMode GeneCounts : output a table of read counts per gene\n", "* --outSAMtype BAM Unsorted : output an unsort BAM file\n", "* --outSAMunmapped Within : included unmapped reads in the BAM file\n", "* --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.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### STAR Output\n", "So what happened? Let's take a look . . ." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam\n", "21_2019_P_M1_S21_L001_R1_short_introns_Aligned.sortedByCoord.out.bam.bai\n", "21_2019_P_M1_S21_L001_R1_short_introns_Log.final.out\n", "21_2019_P_M1_S21_L001_R1_short_introns_Log.out\n", "21_2019_P_M1_S21_L001_R1_short_introns_Log.progress.out\n", "21_2019_P_M1_S21_L001_R1_short_introns_ReadsPerGene.out.tab\n", "21_2019_P_M1_S21_L001_R1_short_introns_SJ.out.tab\n", "21_2019_P_M1_S21_L002_R1_Aligned.out.bam\n", "21_2019_P_M1_S21_L002_R1_Log.final.out\n", "21_2019_P_M1_S21_L002_R1_Log.out\n", "21_2019_P_M1_S21_L002_R1_Log.progress.out\n", "21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab\n", "21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam\n", "21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai\n", "21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out\n", "21_2019_P_M1_S21_L002_R1_short_introns_Log.out\n", "21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out\n", "21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab\n", "21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab\n", "21_2019_P_M1_S21_L002_R1_SJ.out.tab\n", "genome_Log.out\n" ] } ], "source": [ "ls ${STAR_OUT}" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Aligned.out.bam <==\n", "\u001f\u0004\u0000\u0000\u0000\u0000\u0000�\u0006\u0000BC\u0002\u0000�\u0002��Mo�@\u0010���C\u001d", "$�@�H�X\\!���_�S�T�JM�\u0012\u0017�M���^����\n", "�\u0013\u0017\u0004?�\u0013\u0012\u0012\u0012g�\t=�\u0013�r�#�q�R�\n", "�#�M\u0018�AX�qLlӪ\u001aVq,�1gy��^�4�Z{p_\u001b(���Z�2b2V2O\u00024\u0019�\u0010�\t�)�\u001d", "�Lk\u0002�\u001c", "���d&0�.��,`��0������v��u�5\u0019�p��#�ڬb\f", "\u0018qm\u0017�e-�\u0000�Ř=�\n", "6:�Z�W���1?��a5��]�����z�\t/\u001e", "q=�\t�b�\u0007���\u001b�D\u001d", "{:\u00050�\"��r���D\t�q�t?\u0010�^�<1�a\u0012����\n", "�\u00122\u0011\u0005ө������\f", "�f����Y�e\u0012F\u0011\u001f��A*wS�\u0012�b����\u0011�Nk\u0018S��\u0018�F�ۀ���@5��g7��@��q&\u000e�i\u001e", "\u001bB��=�d>V/�x?��ᳳN�� �U���\u0019TƠ���)ח:�B\u000fD\u001a'���\n", "\u0017�\u001bgB�ިO�.��i��\u0010�ł���{Z��D\u001f\u001e", "��\u0013\n", "�5YY�y�\u001b�ϗ��^����U]�j&��s�9������\u0005\u001440=b\u0018�\u0006;\u0018�\u001a�\u001d", "\u0011#�ĮMEKGC\u000b", "QT\u0010;a\\.pF1�AŏN\u0000����cl\u0003\u0013�J�a���Y�0���Uu�{�ǩ�Y+����ԩ�{������~��\u0007g��\u0015~����0{��^\u001d", "������Twp:{�Y�͏�����8|F�����s^b��x�|�����x��������G�fy^�Ɩ���&��8I��\u0013[$u�DM\u0012Em\u0016��ud�q՝������'��������\u0017Vw�T=�Λ�\u0017���W����\u0007\u0019�c���\u000f�`z\u0003�n�\u0019ӒC��\u001d", "�(lӠ��1\u0019�\u0002��aS�eX��d�i�*�L?�S#<��i̙�j���C�\u0015���)tZ�\u000e�{�t�0\u0015�i�,HB�%a��IP\u0015U�f6.m��\"�y\u0012$n0P�H��P:���Q�^��\u0015\u0011�ʕk�r%��]�~���kׯ�\u000fl�ӣ��_x����I�櫽�\t&\u0013DE\u0015\u0007aY�Q�\u0014�^\u001b�E\u001c", "�MY��j��$�q��咇\u0003��lwq�N�\n", "�[���e���O�c��.>�0�\u000eVLJ\u001e", "�;4��06a֤Un�,\f", "�866��&/�0,�\"/��o��&�#�bl��\u0012U����܍ۿ�k����k\u0010\u000b", ":e#D��l�\u0000�\u001c", "?�\u001f<\u0019D�h\n", "�V\u0007\u000e���\u000b", "�&,�6o\"�D��\u0011G�Vi\u001bۤ��\"�aRY/ ����\u001c", "�\u001e", "��{\u0007�\u0004��\u000f8ݣ#HN@�\u001ap\u001e", "���\u0003W`�&�\u0019\u0001�<>�@��\u001e", "2\u0011*\u001d", "\u001d", "��J\u0011���4\f", "3\u001b ��&(�\u0017m�ب��\u0018U)׃\u0012X@�׮�!�h�t�N��\u0015��\n", "\u000f8tW�z1�\u0003�̛\u0012�9��\u0013�W /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.final.out <==\n", " Started job on |\tJun 26 15:10:01\n", " Started mapping on |\tJun 26 15:10:05\n", " Finished on |\tJun 26 15:11:40\n", " Mapping speed, Million of reads per hour |\t92.30\n", "\n", " Number of input reads |\t2435761\n", " Average input read length |\t75\n", " UNIQUE READS:\n", " Uniquely mapped reads number |\t2354695\n", " Uniquely mapped reads % |\t96.67%\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.out <==\n", "STAR version=STAR_2.5.2b\n", "STAR compilation time,server,dir=\n", "##### DEFAULT parameters:\n", "versionSTAR 20201\n", "versionGenome 20101 20200 \n", "parametersFiles - \n", "sysShell -\n", "runMode alignReads\n", "runThreadN 1\n", "runDirPerm User_RWX\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_Log.progress.out <==\n", " Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped\n", " M/hr number length unique length MMrate multi multi+ MM short other\n", "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%\n", "ALL DONE!\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_ReadsPerGene.out.tab <==\n", "N_unmapped\t46060\t46060\t46060\n", "N_multimapping\t35006\t35006\t35006\n", "N_noFeature\t12466\t2145291\t18783\n", "N_ambiguous\t203327\t820\t316\n", "CNAG_04548\t0\t0\t0\n", "CNAG_07303\t0\t0\t0\n", "CNAG_07304\t6\t0\t6\n", "CNAG_00001\t0\t0\t0\n", "CNAG_07305\t1\t0\t1\n", "CNAG_00002\t51\t0\t51\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam <==\n", "�F��T��\\��i����D\u0004i���(\u0019�\u000e�{y�ׇ��L�v33���S�[��f�nJ\u000b", "^B�;�\u001bݝ�ɬ}A�\u0007я�I\u0010\u0004��\u0011z�#�\u0019������\u0016\u000b", "B\u0003a����̓�Fp��u�\u0016�,\u0019\u001b�:��F�V�i�KD,ye!|f�j�Xm�)e�f�`�\u0004h3R\":%ĥ����Հ9�ClR0\u001b\u0018�m\u001f�>\u0007��V8�u�ٔ�n�`\u001e", "�\u001ac\u000e.��|���\u0001#��c�<\u0006�\u000b", ";��v\t�ޡK\u0018;�p\u0014ϳ�KK\u0006G�\u001c", "׳�}p\u0014fS��*\u000b", "�ece�\u001e", "�7:F���\n", "`�E:�KIfZC%�\u0017��XX�i��ʻY,�Ck+I\u0013я\u0012!���\u0005z25��d��\u0015q�j�%�1�E�8�;�EID1aQ;\n", "H\u0014��*�4cR�z�`�\u000e\u000eLxN��L��X���D\u001c", "$�ӱn�I'rji�c��x?ٻ��r\u0019g��>_�0ͤv�<*V\u000b", "\u001b�����b#P(T\u0016�[�o�I?���i?\u0019I���\u0014A���w�:\u001c", "���`h�nn��O}\u0018��$\u0016\u0012�x�\u000b", "�L'B�E��M��<��V*ԯ�\\3&9����g\u001b%��M�˙a9=��咣r~JTF.3\u001e", "�\\ҩ̜������\u0011�\u0000Q�=��Wn!�������\u001f�kJ\u0005\u0012m�\u001aR���}P5�kF+\u001b��8\u000b", "�A?�΁rћowAy��V>���\u0015C�n��\u0004�׏f�$���{ZR�����54��H��N�&�\u000b", "�Ъ�.\u0006\u0000\u0000\u001f\u0004\u0000\u0000\u0000\u0000\u0000�\u0006\u0000BC\u0002\u0000e=�}{�,[VV�{�0�h���:�z�9Ֆڀ@�_�s�ny��\u0019�+�0*��#�\t����11��Ҵ<�����\n", "�\u0006\u0013�\u0000\u0006u�1b4����׮]���Uu�ĩ{�U��O�W��k����{���\u001c", "���F��{�����{Q�\u001b�������b?�\u001fW��1u�{5}�\u0007���X.\u0016��b�����~5��K�Q\u0012\u0004E\u0016dyz\u000e��\u001c", "&�)+�I0\t�s\u0011M��y���o<�\u0011_\u001c", "�E\u000f�>�ǟ���_y)�+/&_��{��\u000f�2�\u0011�h�݀��#\u0003�}]!-\u0017���߮���b�\u001a}\f", "����ir�ϓ�)���i��t��i��8\u001d", "G�N� ��{�\u0002@z; }\u0016 }�\u000f��\u0004�,K�!-\u0000i{\\�W��j7�4ˣl�f�8�&�I�f��l\u001a\u0015i:�r�F�l�N�Cz��n\u000eKΓXY�o�~R=��cI }²�.@R��\u0005i\u0007�V��q��n\u0005���\u0014\u0015���HNI���p�M�I�g�)�\u0014��$I���'-3OP&x\n", "����)�7\t���[��\u0004�?��6���\u0010�0?CA�I2��\u0019�l�AI��4:EA\u0000\u0019\f", "'Y;S�Z�-�B�����c�!j�'\u001a�@-��������j��DM�\u000e\u0007u\u001d", "��\u001d", "OF)~ڢj��S���a�ܬ�#�\u000b", "��Z4�mG�q\u0010L�i���(\t�b�J_�c\u001e", "���MP���-����\u0007���\u0002�q�8,�����Ί�*o��`�\u0007�\u001a��\u001a\u0004i��\u001f���\u001f��\u001a��\u0015v\u000b", "=S��(,\f", "�{�m�w�\u0012�>�7\"�=}��\u0000%�\u0004�5�N��f�[\u001b�|�nb̻\t̻S8���\u0018\n", "#���\u0015��ڕ��\u0012�a=z��\u0012\u0000�m�~t\u0010(XF\u0007\u0018\u0010�\u0005T�a=7�����0�a.����\u0018%���f��\u0006%2؎\u0000=u㠒x��О`\u0015-�[���0L=\u000fP�j��u� w������Z-�����򸜯�\u0006�5n>m���&�����G\"Ӧ^�ug+��\u001c", "�|�_ÉZ�a\"��9}�8\u000e�I\u0012'ERde\u0010���FPx콒4+2|�K�-�>=��\u0001�\u001d", "hW\u000f�NԿ\u0001�\u0001*��o=��n��\u0007\u0006�E\u0018O\u001e", "\u0011�W��T�\n", "�6�)\f", "�����u�<[����c\u000b", "�e\u000e� `\t�������\u000f}ᓠ���?^A��ӂ·\\��#z`c�O�K�`\u0005\u000e.\u0014����\\F�d*\\U�U�Gk��P\u001e", "#x>��D���OX�GL�:��:�d\f", "�\f", "�����\u0013ކ=PHo\u0004>\t\u000e\u000b", "��Wsw\f", "�ԏ�\u0011����\\\u001e", "���~+\n", "��\u001bđև�R�\u001f}�����\\�R�29\u0004�C�{�g�G���Ef��/�����\u0016�\u001e", "5p�TQ��kS=\u0015\u0005��55�b�[\u0018�ׅ� ��2�z�II \u0000?l*\n", "�{��!�\u000e���Zn\u0016;t‹�A�DC1�S\u001a�e\u001e", "��s\u0018��� x��\"=<|�K�z�|C��\u0006\tP��#�{�F\u0004}��ٜ�#<\u000eL��Q�g\u0010�\u0017�\t�\u0012�v�v{6���O�a:�\u001d", "��br��\u001e", "C*��9��\u0017������7yz�g@���s�\u0018\u0016��at\u001f]wt�\f", "���\f", "��\f", "�fE\u0014���>UHlS�\u000f�+�[@z��$��ǭ���S��_\u0018��CO@�\u0018i��I\\����\u0007\u001e", "�F�����,I�I��)>'Y\u001e", "N�<���$/\u0011�@\u00183J��S�.�\tP9ƾx\u0001\u0001ca�C\u001a&`џb/������߻\u0000˨����\u0013�8(��r�H�|]�~I���8,ð\f", "���9��E8��0��cj��z\u0004�c4\u0013��l\u0019�#�\u0017��Ǐ�t��-XF��3Ջ�\u0015���|�;\u001e", "�\u001aqn�v\n", "B#�x\u0001�\t�:��A|յ����\u001c", "���[,w{c�u���\u001d", "y*��\t4�T�U��I��\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Aligned.sortedByCoord.out.bam.bai <==\n", "�J\u0015\u0000\u0000\u0000M\u0002\u0000\u0000\u0004\u0000\u0000\u0000����(\u0000\u0000\u0000�c�1+\u0000\u0000\u0000���B,\u0000\u0000\u0000�Έ+-\u0000\u0000\u0000\u0011rpT/\u0000\u0000\u0000+�pT/\u0000\u0000\u0000\u001at�a0\u0000\u0000\u0000�rU�0\u0000\u0000\u0000`\u0012\u0000\u0000\u0001\u0000\u0000\u0000�0ϱ\u001e", "\u0000\u0000\u0000��\u0018�\u001f\u0000\u0000\u0000^\u0012\u0000\u0000\u0001\u0000\u0000\u0000��N�\u001a\u0000\u0000\u0000�7^>\u001d", "\u0000\u0000\u0000n\u0012\u0000\u0000\u0001\u0000\u0000\u0000�Έ+-\u0000\u0000\u0000\u0011rpT/\u0000\u0000\u0000o\u0012\u0000\u0000\u0001\u0000\u0000\u0000+�pT/\u0000\u0000\u0000\u001at�a0\u0000\u0000\u0000l\u0012\u0000\u0000\u0001\u0000\u0000\u0000�c�1+\u0000\u0000\u0000���B,\u0000\u0000\u0000r\u0012\u0000\u0000\u0001\u0000\u0000\u0000s\u0007p�1\u0000\u0000\u0000�Y�;5\u0000\u0000\u0000z\u0012\u0000\u0000\u0001\u0000\u0000\u0000���N<\u0000\u0000\u0000q\u0017Mj=\u0000\u0000\u0000~\u0012\u0000\u0000\u0001\u0000\u0000\u0000^�\u0001LC\u0000\u0000\u0000Z\u001c", "�`D\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�j�\u0007E\u0000\u0000\u0000O\u0018[�F\u0000\u0000\u0000}\u0012\u0000\u0000\u0001\u0000\u0000\u0000f �]>\u0000\u0000\u0000�c\u0001LC\u0000\u0000\u0000N\u0002\u0000\u0000\u0002\u0000\u0000\u00004yU�0\u0000\u0000\u0000s\u0007p�1\u0000\u0000\u0000\u001f��M6\u0000\u0000\u0000��av9\u0000\u0000\u0000O\u0002\u0000\u0000\u0004\u0000\u0000\u00007���;\u0000\u0000\u0000���N<\u0000\u0000\u0000q\u0017Mj=\u0000\u0000\u0000f �]>\u0000\u0000\u0000X,\u0001LC\u0000\u0000\u0000^�\u0001LC\u0000\u0000\u0000Z\u001c", "�`D\u0000\u0000\u0000�j�\u0007E\u0000\u0000\u0000x\u0012\u0000\u0000\u0001\u0000\u0000\u0000��av9\u0000\u0000\u0000[���;\u0000\u0000\u0000s\u0012\u0000\u0000\u0001\u0000\u0000\u0000�Y�;5\u0000\u0000\u0000\u001f��M6\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000nE_pG\u0000\u0000\u0000�ڣ�I\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000!�jvM\u0000\u0000\u0000\u001f��JQ\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u001f��JQ\u0000\u0000\u0000�.K�T\u0000\u0000\u0000P\u0002\u0000\u0000\u0002\u0000\u0000\u0000�2[�F\u0000\u0000\u0000nE_pG\u0000\u0000\u0000�ڣ�I\u0000\u0000\u0000�(L�L\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u0010�|KV\u0000\u0000\u0000�~:�W\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�:�W\u0000\u0000\u0000Vq��X\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000Vq��X\u0000\u0000\u0000���Z\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000���Z\u0000\u0000\u0000k.h�^\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000>�T�k\u0000\u0000\u0000�Y��l\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000W�\n", "\\`\u0000\u0000\u0000\u0013\u0011+xc\u0000\u0000\u0000Q\u0002\u0000\u0000\u0003\u0000\u0000\u0000�*L�L\u0000\u0000\u0000!�jvM\u0000\u0000\u0000�.K�T\u0000\u0000\u0000\u0010�|KV\u0000\u0000\u0000�~:�W\u0000\u0000\u0000�:�W\u0000\u0000\u0000R\u0002\u0000\u0000\u0002\u0000\u0000\u0000I���_\u0000\u0000\u0000W�\n", "\\`\u0000\u0000\u0000I\u000e'�e\u0000\u0000\u0000��բf\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000k.h�^\u0000\u0000\u0000I���_\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000��բf\u0000\u0000\u0000��g\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u0013\u0011+xc\u0000\u0000\u0000I\u000e'�e\u0000\u0000\u0000J\u0000\u0000\u0000\u0005\u0000\u0000\u0000���Z\u0000\u0000\u0000���Z\u0000\u0000\u0000��g\u0000\u0000\u0000�\u000b", "��g\u0000\u0000\u0000<� \u0010o\u0000\u0000\u0000�� \u0010o\u0000\u0000\u0000󉱎x\u0000\u0000\u0000����x\u0000\u0000\u000050�d�\u0000\u0000\u000005�d�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000���\u0000h\u0000\u0000\u0000\n", "�\u0001\u001e", "i\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000^�\u0001\u001e", "i\u0000\u0000\u0000��T�k\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000<&xPq\u0000\u0000\u0000Jʹ\\r\u0000\u0000\u0000S\u0002\u0000\u0000\u0004\u0000\u0000\u0000�\u000b", "��g\u0000\u0000\u0000���\u0000h\u0000\u0000\u0000\n", "�5�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000{y\u0019l~\u0000\u0000\u00005T�̂\u0000\u0000\u0000U\u0002\u0000\u0000\u0004\u0000\u0000\u0000{�nz\u0000\u0000\u0000�m6Q{\u0000\u0000\u0000��2`|\u0000\u0000\u0000馷\u0016}\u0000\u0000\u0000Zp\u0019l~\u0000\u0000\u0000{y\u0019l~\u0000\u0000\u00005T�̂\u0000\u0000\u000050�d�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�m6Q{\u0000\u0000\u0000��2`|\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�\u000e�5�\u0000\u0000\u0000,lYÎ\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000׉�\u0000\u0000\u0000�@ҋ�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000q�t\u0007�\u0000\u0000\u0000�B�\u0019�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�B�\u0019�\u0000\u0000\u0000\u0000\u0001\u0000\u0000\u000005�d�\u0000\u0000\u0000�0>\u001f�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�0>\u001f�\u0000\u0000\u0000�\n", "��?�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\n", "�J\u0015\u0000\u0000\u0000΢�J\u0015\u0000\u0000\u0000S�j�\u0016\u0000\u0000\u0000�\u0002R�\u0017\u0000\u0000\u0000f\"N�\u001a\u0000\u0000\u0000\u0000�N�\u001a\u0000\u0000\u0000�7^>\u001d", "\u0000\u0000\u0000����\u001e", "\u0000\u0000\u0000��\u0018�\u001f\u0000\u0000\u0000 \u001c", "�� \u0000\u0000\u00001�\u0010�!\u0000\u0000\u0000�\tx�&\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000 �\u0006J�\u0000\u0000\u0000\u0004q:1�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u0003\u0018�ޚ\u0000\u0000\u0000����\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u0006�ƹ�\u0000\u0000\u0000d�d��\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000=x���\u0000\u0000\u0000`\u0017�ի\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�:1�\u0000\u0000\u0000G�g;�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000\u0013�g;�\u0000\u0000\u0000�=ƹ�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�,�(�\u0000\u0000\u0000JAD��\u0000\u0000\u0000X\u0002\u0000\u0000\u0005\u0000\u0000\u0000��xٜ\u0000\u0000\u0000 �\u0006J�\u0000\u0000\u0000\u0004q:1�\u0000\u0000\u0000�:1�\u0000\u0000\u0000G�g;�\u0000\u0000\u0000\u0013�g;�\u0000\u0000\u0000�6ƹ�\u0000\u0000\u0000\u0006�ƹ�\u0000\u0000\u0000d�d��\u0000\u0000\u0000^e�$�\u0000\u0000\u0000Y\u0002\u0000\u0000\u0004\u0000\u0000\u0000^e�$�\u0000\u0000\u0000=x���\u0000\u0000\u0000`\u0017�ի\u0000\u0000\u00008\u0018�ի\u0000\u0000\u0000��&T�\u0000\u0000\u0000�,�(�\u0000\u0000\u0000JAD��\u0000\u0000\u0000�d�Y�\u0000\u0000\u0000Z\u0002\u0000\u0000\u0002\u0000\u0000\u0000�d�Y�\u0000\u0000\u0000�n���\u0000\u0000\u0000�����\u0000\u0000\u0000�ζs�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u00008\u0018�ի\u0000\u0000\u0000��&T�\u0000\u0000\u0000�\u0012\u0000\u0000\u0001\u0000\u0000\u0000�n���\u0000\u0000\u0000�����\u0000\u0000\u0000J�\u0000\u0000\u0002\u0000\u0000\u0000\u0000\u0000�\u0002\u0000\u0000\u0000\u0000�ζs�\u0000\u0000\u0000!�\u0004\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000\u0000�\u0000\u0000\u0000\u0000\u0000�\u0002\u0000\u0000\u0000\u0000�\u001b\u000b", "@\u0000\u0000\u0000\u0000\u0017�%�\u0000\u0000\u0000\u0000�r�h\u0001\u0000\u0000\u0000���\u0001\u0000\u0000\u0000����\u0003\u0000\u0000\u0000����\u0003\u0000\u0000\u0000y\u0012�/\u0006\u0000\u0000\u0000�w��\u0006\u0000\u0000\u0000��#W\u0007\u0000\u0000\u0000��k�\t\u0000\u0000\u0000:+�j\u000b", "\u0000\u0000\u0000��դ\u000b", "\u0000\u0000\u0000�\u0000E�\n", "�6'\u0000\u0000\u0000�y\u000fz'\u0000\u0000\u0000��8�'\u0000\u0000\u0000���q(\u0000\u0000\u0000����(\u0000\u0000\u0000�@��)\u0000\u0000\u00002��*\u0000\u0000\u00004b�1+\u0000\u0000\u0000���B,\u0000\u0000\u0000�Έ+-\u0000\u0000\u0000\u0011rpT/\u0000\u0000\u0000\u001at�a0\u0000\u0000\u0000�rU�0\u0000\u0000\u0000�\u0004p�1\u0000\u0000\u0000�Y�;5\u0000\u0000\u0000\u001f��M6\u0000\u0000\u0000�\u0010��6\u0000\u0000\u0000��\u0004�7\u0000\u0000\u0000XE��8\u0000\u0000\u0000��av9\u0000\u0000\u0000[���;\u0000\u0000\u0000���N<\u0000\u0000\u0000q\u0017Mj=\u0000\u0000\u0000����=\u0000\u0000\u0000�\u001f�]>\u0000\u0000\u0000X,\u0001LC\u0000\u0000\u0000Z\u001c", "�`D\u0000\u0000\u0000Ff�\u0007E\u0000\u0000\u0000O\u0018[�F\u0000\u0000\u0000nE_pG\u0000\u0000\u0000�ڣ�I\u0000\u0000\u0000��9'J\u0000\u0000\u0000�ğjJ\u0000\u0000\u0000\u001e", "�*�J\u0000\u0000\u0000s\u0014�`K\u0000\u0000\u0000��w�K\u0000\u0000\u0000�(L�L\u0000\u0000\u0000!�jvM\u0000\u0000\u0000\u001f��JQ\u0000\u0000\u0000�.K�T\u0000\u0000\u0000L\u000f��U\u0000\u0000\u0000\u0010�|KV\u0000\u0000\u0000�~:�W\u0000\u0000\u0000Vq��X\u0000\u0000\u0000���Z\u0000\u0000\u0000k.h�^\u0000\u0000\u0000I���_\u0000\u0000\u0000\f", "�\n", "\\`\u0000\u0000\u0000\u0013\u0011+xc\u0000\u0000\u0000I\u000e'�e\u0000\u0000\u0000vR�+f\u0000\u0000\u0000��բf\u0000\u0000\u0000��g\u0000\u0000\u0000���\u0000h\u0000\u0000\u0000\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "�5�\u0000\u0000\u0000,lYÎ\u0000\u0000\u0000�5�\u0002�\u0000\u0000\u0000���\u0000\u0000\u0000�@ҋ�\u0000\u0000\u0000`�t\u0007�\u0000\u0000\u0000�B�\u0019�\u0000\u0000\u0000o\u0000\u0000\u0000���gp\u0000\u0000\u0000Z%xPq\u0000\u0000\u0000Jʹ\\r\u0000\u0000\u0000 ���t\u0000\u0000\u0000�\u0015S�v\u0000\u0000\u0000�\u0014�*w\u0000\u0000\u0000y�7\u0011x\u0000\u0000\u0000󉱎x\u0000\u0000\u0000{�nz\u0000\u0000\u0000�m6Q{\u0000\u0000\u0000��2`|\u0000\u0000\u0000ɘ�\u0016}\u0000\u0000\u0000Zp\u0019l~\u0000\u0000\u00005T�̂\u0000\u0000\u0000�\u001f�\u0000\u0000\u0000�\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.final.out <==\n", " Started job on |\tJun 24 09:50:39\n", " Started mapping on |\tJun 24 09:50:39\n", " Finished on |\tJun 24 09:51:17\n", " Mapping speed, Million of reads per hour |\t230.76\n", "\n", " Number of input reads |\t2435761\n", " Average input read length |\t75\n", " UNIQUE READS:\n", " Uniquely mapped reads number |\t2355964\n", " Uniquely mapped reads % |\t96.72%\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.out <==\n", "STAR version=STAR_2.5.2b\n", "STAR compilation time,server,dir=\n", "##### DEFAULT parameters:\n", "versionSTAR 20201\n", "versionGenome 20101 20200 \n", "parametersFiles - \n", "sysShell -\n", "runMode alignReads\n", "runThreadN 1\n", "runDirPerm User_RWX\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_Log.progress.out <==\n", " Time Speed Read Read Mapped Mapped Mapped Mapped Unmapped Unmapped Unmapped Unmapped\n", " M/hr number length unique length MMrate multi multi+ MM short other\n", "ALL DONE!\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_ReadsPerGene.out.tab <==\n", "N_unmapped\t46059\t46059\t46059\n", "N_multimapping\t33738\t33738\t33738\n", "N_noFeature\t12438\t2146527\t18736\n", "N_ambiguous\t203237\t786\t143\n", "CNAG_04548\t0\t0\t0\n", "CNAG_07303\t0\t0\t0\n", "CNAG_07304\t6\t0\t6\n", "CNAG_00001\t0\t0\t0\n", "CNAG_07305\t1\t0\t1\n", "CNAG_00002\t51\t0\t51\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_short_introns_SJ.out.tab <==\n", "1\t9247\t9310\t2\t2\t1\t6\t4\t35\n", "1\t10455\t10506\t1\t3\t0\t7\t8\t24\n", "1\t10609\t10694\t1\t1\t0\t0\t2\t21\n", "1\t11574\t11622\t2\t2\t0\t0\t1\t21\n", "1\t11574\t11637\t2\t2\t0\t0\t2\t32\n", "1\t12904\t12991\t1\t1\t1\t4\t0\t37\n", "1\t13077\t13136\t1\t1\t1\t10\t0\t37\n", "1\t13256\t13310\t1\t1\t1\t7\t0\t37\n", "1\t13480\t13530\t1\t1\t1\t5\t0\t33\n", "1\t13480\t13534\t1\t1\t0\t1\t0\t28\n", "\n", "==> /home/jovyan/work/scratch/bioinf_intro/star_out/21_2019_P_M1_S21_L002_R1_SJ.out.tab <==\n", "1\t9247\t9310\t2\t2\t1\t6\t4\t35\n", "1\t10455\t10506\t1\t3\t0\t7\t7\t24\n", "1\t10609\t10694\t1\t1\t0\t0\t2\t21\n", "1\t11574\t11622\t2\t2\t0\t0\t1\t21\n", "1\t11574\t11637\t2\t2\t0\t0\t2\t32\n", "1\t12904\t12991\t1\t1\t1\t4\t0\t37\n", "1\t13077\t13136\t1\t1\t1\t10\t0\t37\n", "1\t13256\t13310\t1\t1\t1\t7\t0\t37\n", "1\t13480\t13530\t1\t1\t1\t5\t0\t33\n", "1\t13480\t13534\t1\t1\t0\t1\t0\t28\n" ] } ], "source": [ "head ${STAR_OUT}/21_2019_P_M1_S21_L002_R1*" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "STAR generates several files for each FASTQ:\n", "* Log.out : lots of details of the run, including all parameters used\n", "* Log.final.out : Important summary statistics\n", "* ReadsPerGene.out.tab : Count table, the main thing we are interested in\n", "* SJ.out.tab : All splice junctions, including ones from the GTF and novel junctions discovered by STAR\n", "* Log.progress.out: run statistics updated during run, not so interesting at the end\n", "\n", "Let's take a closer look at Log.final.out" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Started job on |\tJun 26 15:10:01\n", " Started mapping on |\tJun 26 15:10:05\n", " Finished on |\tJun 26 15:11:40\n", " Mapping speed, Million of reads per hour |\t92.30\n", "\n", " Number of input reads |\t2435761\n", " Average input read length |\t75\n", " UNIQUE READS:\n", " Uniquely mapped reads number |\t2354695\n", " Uniquely mapped reads % |\t96.67%\n", " Average mapped length |\t75.36\n", " Number of splices: Total |\t610596\n", " Number of splices: Annotated (sjdb) |\t596312\n", " Number of splices: GT/AG |\t603894\n", " Number of splices: GC/AG |\t5846\n", " Number of splices: AT/AC |\t112\n", " Number of splices: Non-canonical |\t744\n", " Mismatch rate per base, % |\t0.16%\n", " Deletion rate per base |\t0.01%\n", " Deletion average length |\t1.35\n", " Insertion rate per base |\t0.00%\n", " Insertion average length |\t1.11\n", " MULTI-MAPPING READS:\n", " Number of reads mapped to multiple loci |\t35006\n", " % of reads mapped to multiple loci |\t1.44%\n", " Number of reads mapped to too many loci |\t1263\n", " % of reads mapped to too many loci |\t0.05%\n", " UNMAPPED READS:\n", " % of reads unmapped: too many mismatches |\t0.00%\n", " % of reads unmapped: too short |\t1.83%\n", " % of reads unmapped: other |\t0.00%\n", " CHIMERIC READS:\n", " Number of chimeric reads |\t0\n", " % of chimeric reads |\t0.00%\n" ] } ], "source": [ "cat ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_Log.final.out" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Sanity Check: Number of Reads\n", "Are the number of input reads what we expect? Let's look at how many reads are in the input FASTQ" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "2435761\n" ] } ], "source": [ "zcat $TRIMMED/21_2019_P_M1_S21_L002_R1_001.trim.fastq.gz | awk '{s++}END{print s/4}' " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Sanity Check: Unmapped Reads\n", "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.\n", "\n", "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](https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ ">NB501800:327:HF27FBGXB:2:11202:18025:8102\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:20718:8137\n", "TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTGTAAAAGTGGGGGTTGTTTTTTATTTTTTTGTAGATTTAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:9226:8227\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:17164:8437\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:25989:8488\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTAACTCGTATGCCGTCTTATGCTTGAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:11951:8497\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:11505:8555\n", "ATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:21516:8563\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAATAAG\n", ">NB501800:327:HF27FBGXB:2:11202:8434:8590\n", "GATCGGAAGAGCACACGTCTGAACTCCAGTCACCGCTCATTATCTCGTATGCCGTCTTCTGCTTGAAAAAAAA\n", ">NB501800:327:HF27FBGXB:2:11202:12081:8652\n", "GATCGAAAGACCACACGTCTGAACACCAGTCACAGCTCAATATCACGTATGAAGTATTATGCATTAAAAAAAAA\n", "[bam2fq_mainloop] Error writing to FASTx files.: Broken pipe\n", "[M::bam2fq_mainloop] discarded 0 singletons\n", "[M::bam2fq_mainloop] processed 71 reads\n", "samtools bam2fq: Error closing STDOUT: Broken pipe\n" ] } ], "source": [ "samtools fasta -f 4 ${STAR_OUT}/21_2019_P_M1_S21_L002_R1_Aligned.out.bam | head -n20" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### MultiQC\n", "MultiQC also works with STAR reports, so let's try it!" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/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.\n", " configs = yaml.load(f)\n", "/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.\n", " sp = yaml.load(f)\n", "[INFO ] multiqc : This is MultiQC v1.7\n", "[INFO ] multiqc : Template : default\n", "[INFO ] multiqc : Searching '/home/jovyan/work/scratch/bioinf_intro/star_out'\n", "\u001b[?25lSearching 21 files.. [####################################] 100% \u001b[?25h\n", "[INFO ] star : Found 3 reports and 3 gene count files\n", "[INFO ] multiqc : Compressing plot data\n", "[INFO ] multiqc : Report : ../../scratch/bioinf_intro/star_out/multiqc_report.html\n", "[INFO ] multiqc : Data : ../../scratch/bioinf_intro/star_out/multiqc_data\n", "[INFO ] multiqc : MultiQC complete\n" ] } ], "source": [ "multiqc ${STAR_OUT} --outdir ${STAR_OUT}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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 :" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/home/jovyan/work/scratch/bioinf_intro/star_out\n" ] } ], "source": [ "echo ${STAR_OUT}" ] } ], "metadata": { "kernelspec": { "display_name": "Bash", "language": "bash", "name": "bash" }, "language_info": { "codemirror_mode": "shell", "file_extension": ".sh", "mimetype": "text/x-sh", "name": "bash" } }, "nbformat": 4, "nbformat_minor": 1 }