text:
MrBayes runs as single-core, but hangs when multi-core
I compiled MrBayes v3.2.7 from github, with the configuration options ./configure --without-beagle --enable-mpi. The program is called mb.
I tried one of the included example files "sceloporus.nex". I find that if I run mb sceloporus.nex, then it appears to make progress, and writes files to disk.
I used MPI with 1 core, and it behaves as before, making progress, writing files to disk.
But if I run with more than 1 core, using MPI, with mpirun -np 2 mb sceloporus.nex, then it only makes it through some of the analysis setup, doesn't write any files to disk, and just hangs. If I look at top, I see activity on the number of cores I specified. I finally had to manually terminate the program.
Is there something I should be doing differently to have MrBayes properly run in parallel?
Here is a copy of my terminal from these 3 attempts:
[ptitle@iu mbTest]$ ls
sceloporus.nex
[ptitle@iu mbTest]$ ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (1 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
      Taxon  11 -> AZlaDGM856
      Taxon  12 -> AZmaDGM984
      Taxon  13 -> AZmaDGM992
      Taxon  14 -> AZmaRM4401
      Taxon  15 -> AZmaRM4402
      Taxon  16 -> AZmaRM4403
      Taxon  17 -> AZmoGM1045
      Taxon  18 -> AZmoDGM492
      Taxon  19 -> AZmoDGM699
      Taxon  20 -> AZmoDGM704
      Taxon  21 -> AZmoRM4394
      Taxon  22 -> AZpiKWS238
      Taxon  23 -> AZyuDGM826
      Taxon  24 -> AZyuDGM838
      Taxon  25 -> AZyuDGM852
      Taxon  26 -> AZyuDGM858
      Taxon  27 -> AZyuDGM865
      Taxon  28 -> CAfrP25750
      Taxon  29 -> CAimDGM534
      Taxon  30 -> CAimDGM558
      Taxon  31 -> CAimDGM694
      Taxon  32 -> CAimKWS170
      Taxon  33 -> CAinCRF530
      Taxon  34 -> CAinoDGM74
      Taxon  35 -> CAinCRF529
      Taxon  36 -> CAkeDS2SP7
      Taxon  37 -> CAkeP26556
      Taxon  38 -> CAla78323
      Taxon  39 -> CAlaDAW220
      Taxon  40 -> CAlaM14610
      Taxon  41 -> CAlaM23289
      Taxon  42 -> CAmoJQR132
      Taxon  43 -> CArvJOS138
      Taxon  44 -> CArvESA441
      Taxon  45 -> CAsaBUR167
      Taxon  46 -> CAsaMCFMAG
      Taxon  47 -> CAsarnMCC
      Taxon  48 -> CAsa150079
      Taxon  49 -> CAsa150088
      Taxon  50 -> CAsa179886
      Taxon  51 -> CAsaM19838
      Taxon  52 -> CAsdDGM691
      Taxon  53 -> CAsdDGM748
      Taxon  54 -> CAsdSWT701
      Taxon  55 -> MXBCDGM474
      Taxon  56 -> MXBCDGM475
      Taxon  57 -> MXBCDGM486
      Taxon  58 -> MXchRM4860
      Taxon  59 -> MXduA83134
      Taxon  60 -> MXsoM14488
      Taxon  61 -> MXsoP26449
      Taxon  62 -> MXso210345
      Taxon  63 -> NMdoTWR230
      Taxon  64 -> NMedDGM938
      Taxon  65 -> NMhiU48819
      Taxon  66 -> NMluGM947
      Taxon  67 -> NMluDGM948
      Taxon  68 -> NMsiDGM891
      Taxon  69 -> NMsoDGM900
      Taxon  70 -> NMsoDGM904
      Taxon  71 -> NMsoDGM905
      Taxon  72 -> NMvaDGM924
      Taxon  73 -> NVchDGM612
      Taxon  74 -> NVch162077
      Taxon  75 -> NVclDGM296
      Taxon  76 -> NVclDGM678
      Taxon  77 -> NVesDGM210
      Taxon  78 -> NVlyDGM610
      Taxon  79 -> NVwaDGM608
      Taxon  80 -> UTemGM1058
      Taxon  81 -> UTgrDGM156
      Taxon  82 -> UTgrDGM169
      Taxon  83 -> UTgrDGM175
      Taxon  84 -> UTkaGM1034
      Taxon  85 -> UTsnDM142
      Taxon  86 -> UTsnDG043
      Taxon  87 -> UTsnnDGM645
      Taxon  88 -> UTwsGM261
      Taxon  89 -> UTwsGM263
      Taxon  90 -> UTwsnDGM751
      Taxon  91 -> UTwsJRM4437
      Taxon  92 -> hunsiM376
      Taxon  93 -> hunsiM377
      Taxon  94 -> hunsiRWM619
      Taxon  95 -> lickGM378
      Taxon  96 -> lineatulus
      Taxon  97 -> orcuOM14377
      Taxon  98 -> orcuOM14453
      Taxon  99 -> orcuOM14673
      Taxon 100 -> orcuRWM798
      Taxon 101 -> orcuS201108
      Taxon 102 -> orcuS201124
      Taxon 103 -> zostZ16290
      Taxon 104 -> zosV161292
      Taxon 105 -> zosV161293
      Taxon 106 -> zosV161311
      Taxon 107 -> zosWM1890
      Taxon 108 -> zosGM364
      Taxon 109 -> zosGM365
      Taxon 110 -> zosGM387
      Taxon 111 -> zosGM393
      Taxon 112 -> zosRO332
      Taxon 113 -> zosOM37006
      Taxon 114 -> zosWM1601
      Taxon 115 -> clarkii
      Taxon 116 -> edwardti80
      Taxon 117 -> smaragdnus
      Taxon 118 -> graciosus
      Taxon 119 -> grammicus
      Taxon 120 -> maculosus
      Taxon 121 -> bicanthlis
      Taxon 122 -> torquatus
      Taxon 123 -> variabilis
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1544128156
      Setting output file names to "sceloporus.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Excluding character(s)
      Setting Nst to 6
      Setting Rates to Gamma
      Successfully set likelihood model parameters
      Setting Brlenspr to Unconstrained:GammaDir(1.00,1.0000,0.70,0.50)
      Successfully set prior model parameters
      Setting number of generations to 5000000
      Running Markov chain
      MCMC stamp = 9284864458
      Seed = 753696747
      Swapseed = 1544128156
      Model settings:

         Data not partitioned --
            Datatype  = DNA
            Nucmodel  = 4by4
            Nst       = 6
                        Substitution rates, expressed as proportions
                        of the rate sum, have a Dirichlet prior
                        (1.00,1.00,1.00,1.00,1.00,1.00)
            Covarion  = No
            # States  = 4
                        State frequencies have a Dirichlet prior
                        (1.00,1.00,1.00,1.00)
            Rates     = Gamma
                        The distribution is approximated using 4 categories.
                        Shape parameter is exponentially
                        distributed with parameter (1.00).

      Active parameters:

         Parameters
         ---------------------
         Revmat              1
         Statefreq           2
         Shape               3
         Ratemultiplier      4
         Topology            5
         Brlens              6
         ---------------------

         1 --  Parameter  = Revmat
               Type       = Rates of reversible rate matrix
               Prior      = Dirichlet(1.00,1.00,1.00,1.00,1.00,1.00)

         2 --  Parameter  = Pi
               Type       = Stationary state frequencies
               Prior      = Dirichlet

         3 --  Parameter  = Alpha
               Type       = Shape of scaled gamma distribution of site rates
               Prior      = Exponential(1.00)

         4 --  Parameter  = Ratemultiplier
               Type       = Partition-specific rate multiplier
               Prior      = Fixed(1.0)

         5 --  Parameter  = Tau
               Type       = Topology
               Prior      = All topologies equally probable a priori
               Subparam.  = V

         6 --  Parameter  = V
               Type       = Branch lengths
               Prior      = Unconstrained:GammaDir(1.0,1.0000,0.7,0.5)

      Number of chains per MPI processor = 8

      The MCMC sampler will use the following moves:
         With prob.  Chain will use move
            0.94 %   Dirichlet(Revmat)
            0.94 %   Slider(Revmat)
            0.94 %   Dirichlet(Pi)
            0.94 %   Slider(Pi)
            1.89 %   Multiplier(Alpha)
            9.43 %   ExtSPR(Tau,V)
            9.43 %   ExtTBR(Tau,V)
            9.43 %   NNI(Tau,V)
            9.43 %   ParsSPR(Tau,V)
           37.74 %   Multiplier(V)
           13.21 %   Nodeslider(V)
            5.66 %   TLMultiplier(V)

      Division 1 has 624 unique site patterns
      Initializing conditional likelihoods
      Using standard FMA likelihood calculator for division 1 (single-precision)

      Initial log likelihoods and log prior probs for run 1:
         Chain 1 -- -38239.042144 -- 53.554896
         Chain 2 -- -37743.352685 -- 53.554896
         Chain 3 -- -38049.750737 -- 53.554896
         Chain 4 -- -38759.582835 -- 53.554896

      Initial log likelihoods and log prior probs for run 2:
         Chain 1 -- -38348.591921 -- 53.554896
         Chain 2 -- -39602.594468 -- 53.554896
         Chain 3 -- -38874.065643 -- 53.554896
         Chain 4 -- -39227.723762 -- 53.554896

      Using a relative burnin of 25.0 % for diagnostics

      Chain results (5000000 generations requested):

         0 -- [-38239.042] (-37743.353) (-38049.751) (-38759.583) * [-38348.592] (-39602.594) (-38874.066) (-39227.724) (...0 remote chains...)
^C
[ptitle@iu mbTest]$ ls
sceloporus.nex  sceloporus.nex.mcmc  sceloporus.nex.run1.p  sceloporus.nex.run1.t  sceloporus.nex.run2.p  sceloporus.nex.run2.t
[ptitle@iu mbTest]$ mpirun -np 1 ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (1 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
      Taxon  11 -> AZlaDGM856
      Taxon  12 -> AZmaDGM984
      Taxon  13 -> AZmaDGM992
      Taxon  14 -> AZmaRM4401
      Taxon  15 -> AZmaRM4402
      Taxon  16 -> AZmaRM4403
      Taxon  17 -> AZmoGM1045
      Taxon  18 -> AZmoDGM492
      Taxon  19 -> AZmoDGM699
      Taxon  20 -> AZmoDGM704
      Taxon  21 -> AZmoRM4394
      Taxon  22 -> AZpiKWS238
      Taxon  23 -> AZyuDGM826
      Taxon  24 -> AZyuDGM838
      Taxon  25 -> AZyuDGM852
      Taxon  26 -> AZyuDGM858
      Taxon  27 -> AZyuDGM865
      Taxon  28 -> CAfrP25750
      Taxon  29 -> CAimDGM534
      Taxon  30 -> CAimDGM558
      Taxon  31 -> CAimDGM694
      Taxon  32 -> CAimKWS170
      Taxon  33 -> CAinCRF530
      Taxon  34 -> CAinoDGM74
      Taxon  35 -> CAinCRF529
      Taxon  36 -> CAkeDS2SP7
      Taxon  37 -> CAkeP26556
      Taxon  38 -> CAla78323
      Taxon  39 -> CAlaDAW220
      Taxon  40 -> CAlaM14610
      Taxon  41 -> CAlaM23289
      Taxon  42 -> CAmoJQR132
      Taxon  43 -> CArvJOS138
      Taxon  44 -> CArvESA441
      Taxon  45 -> CAsaBUR167
      Taxon  46 -> CAsaMCFMAG
      Taxon  47 -> CAsarnMCC
      Taxon  48 -> CAsa150079
      Taxon  49 -> CAsa150088
      Taxon  50 -> CAsa179886
      Taxon  51 -> CAsaM19838
      Taxon  52 -> CAsdDGM691
      Taxon  53 -> CAsdDGM748
      Taxon  54 -> CAsdSWT701
      Taxon  55 -> MXBCDGM474
      Taxon  56 -> MXBCDGM475
      Taxon  57 -> MXBCDGM486
      Taxon  58 -> MXchRM4860
      Taxon  59 -> MXduA83134
      Taxon  60 -> MXsoM14488
      Taxon  61 -> MXsoP26449
      Taxon  62 -> MXso210345
      Taxon  63 -> NMdoTWR230
      Taxon  64 -> NMedDGM938
      Taxon  65 -> NMhiU48819
      Taxon  66 -> NMluGM947
      Taxon  67 -> NMluDGM948
      Taxon  68 -> NMsiDGM891
      Taxon  69 -> NMsoDGM900
      Taxon  70 -> NMsoDGM904
      Taxon  71 -> NMsoDGM905
      Taxon  72 -> NMvaDGM924
      Taxon  73 -> NVchDGM612
      Taxon  74 -> NVch162077
      Taxon  75 -> NVclDGM296
      Taxon  76 -> NVclDGM678
      Taxon  77 -> NVesDGM210
      Taxon  78 -> NVlyDGM610
      Taxon  79 -> NVwaDGM608
      Taxon  80 -> UTemGM1058
      Taxon  81 -> UTgrDGM156
      Taxon  82 -> UTgrDGM169
      Taxon  83 -> UTgrDGM175
      Taxon  84 -> UTkaGM1034
      Taxon  85 -> UTsnDM142
      Taxon  86 -> UTsnDG043
      Taxon  87 -> UTsnnDGM645
      Taxon  88 -> UTwsGM261
      Taxon  89 -> UTwsGM263
      Taxon  90 -> UTwsnDGM751
      Taxon  91 -> UTwsJRM4437
      Taxon  92 -> hunsiM376
      Taxon  93 -> hunsiM377
      Taxon  94 -> hunsiRWM619
      Taxon  95 -> lickGM378
      Taxon  96 -> lineatulus
      Taxon  97 -> orcuOM14377
      Taxon  98 -> orcuOM14453
      Taxon  99 -> orcuOM14673
      Taxon 100 -> orcuRWM798
      Taxon 101 -> orcuS201108
      Taxon 102 -> orcuS201124
      Taxon 103 -> zostZ16290
      Taxon 104 -> zosV161292
      Taxon 105 -> zosV161293
      Taxon 106 -> zosV161311
      Taxon 107 -> zosWM1890
      Taxon 108 -> zosGM364
      Taxon 109 -> zosGM365
      Taxon 110 -> zosGM387
      Taxon 111 -> zosGM393
      Taxon 112 -> zosRO332
      Taxon 113 -> zosOM37006
      Taxon 114 -> zosWM1601
      Taxon 115 -> clarkii
      Taxon 116 -> edwardti80
      Taxon 117 -> smaragdnus
      Taxon 118 -> graciosus
      Taxon 119 -> grammicus
      Taxon 120 -> maculosus
      Taxon 121 -> bicanthlis
      Taxon 122 -> torquatus
      Taxon 123 -> variabilis
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1544128192
      Setting output file names to "sceloporus.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Excluding character(s)
      Setting Nst to 6
      Setting Rates to Gamma
      Successfully set likelihood model parameters
      Setting Brlenspr to Unconstrained:GammaDir(1.00,1.0000,0.70,0.50)
      Successfully set prior model parameters
      Setting number of generations to 5000000
      Running Markov chain
      MCMC stamp = 9909051894
      Seed = 1261289135
      Swapseed = 1544128192
      Model settings:

         Data not partitioned --
            Datatype  = DNA
            Nucmodel  = 4by4
            Nst       = 6
                        Substitution rates, expressed as proportions
                        of the rate sum, have a Dirichlet prior
                        (1.00,1.00,1.00,1.00,1.00,1.00)
            Covarion  = No
            # States  = 4
                        State frequencies have a Dirichlet prior
                        (1.00,1.00,1.00,1.00)
            Rates     = Gamma
                        The distribution is approximated using 4 categories.
                        Shape parameter is exponentially
                        distributed with parameter (1.00).

      Active parameters:

         Parameters
         ---------------------
         Revmat              1
         Statefreq           2
         Shape               3
         Ratemultiplier      4
         Topology            5
         Brlens              6
         ---------------------

         1 --  Parameter  = Revmat
               Type       = Rates of reversible rate matrix
               Prior      = Dirichlet(1.00,1.00,1.00,1.00,1.00,1.00)

         2 --  Parameter  = Pi
               Type       = Stationary state frequencies
               Prior      = Dirichlet

         3 --  Parameter  = Alpha
               Type       = Shape of scaled gamma distribution of site rates
               Prior      = Exponential(1.00)

         4 --  Parameter  = Ratemultiplier
               Type       = Partition-specific rate multiplier
               Prior      = Fixed(1.0)

         5 --  Parameter  = Tau
               Type       = Topology
               Prior      = All topologies equally probable a priori
               Subparam.  = V

         6 --  Parameter  = V
               Type       = Branch lengths
               Prior      = Unconstrained:GammaDir(1.0,1.0000,0.7,0.5)

      Number of chains per MPI processor = 8

      The MCMC sampler will use the following moves:
         With prob.  Chain will use move
            0.94 %   Dirichlet(Revmat)
            0.94 %   Slider(Revmat)
            0.94 %   Dirichlet(Pi)
            0.94 %   Slider(Pi)
            1.89 %   Multiplier(Alpha)
            9.43 %   ExtSPR(Tau,V)
            9.43 %   ExtTBR(Tau,V)
            9.43 %   NNI(Tau,V)
            9.43 %   ParsSPR(Tau,V)
           37.74 %   Multiplier(V)
           13.21 %   Nodeslider(V)
            5.66 %   TLMultiplier(V)

      Division 1 has 624 unique site patterns
      Initializing conditional likelihoods
      Using standard FMA likelihood calculator for division 1 (single-precision)

      Initial log likelihoods and log prior probs for run 1:
         Chain 1 -- -39098.978535 -- 53.554896
         Chain 2 -- -38054.044640 -- 53.554896
         Chain 3 -- -39425.015629 -- 53.554896
         Chain 4 -- -40157.756230 -- 53.554896

      Initial log likelihoods and log prior probs for run 2:
         Chain 1 -- -40084.553712 -- 53.554896
         Chain 2 -- -39002.725314 -- 53.554896
         Chain 3 -- -38719.752582 -- 53.554896
         Chain 4 -- -39156.129424 -- 53.554896

      Using a relative burnin of 25.0 % for diagnostics

      Chain results (5000000 generations requested):

         0 -- [-39098.979] (-38054.045) (-39425.016) (-40157.756) * [-40084.554] (-39002.725) (-38719.753) (-39156.129) (...0 remote chains...)
^C[ptitle@iu mbTest]$ mpirun -np 2 ~/MrBayes/src/mb sceloporus.nex

                            MrBayes pre-3.2.7 x86_64

                      (Bayesian Analysis of Phylogeny)

                             (Parallel version)
                         (2 processors available)

              Distributed under the GNU General Public License

               Type "help" or "help <command>" for information
                     on the commands that are available.

                   Type "about" for authorship and general
                       information about the program.

   Executing file "sceloporus.nex"
   DOS line termination
   Longest line length = 1632
   Parsing file
   Expecting NEXUS formatted file
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 123 taxa and 1606 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> AZYuJAS289
      Taxon   2 -> AZYuJAS290
      Taxon   3 -> AZlaJAS294
      Taxon   4 -> AZcoTBP271
      Taxon   5 -> AZcoGM1054
      Taxon   6 -> AZcoDGM994
      Taxon   7 -> AZgiP26438
      Taxon   8 -> AZgrGM1012
      Taxon   9 -> AZgrGM1013
      Taxon  10 -> AZgrGM1024
^C[ptitle@iu mbTest]$

It is not trivial getting MB to run in parallel because from my experience it required granting privileges to regular users. I used this via MPI (MPICH) on a Slurm/SGE cluster.
Essentially, this is coarse grain parallel processing.
It is certainly worth parallelizing MB because each chain can be run on a different nodes and avoids RAM bottlenecks on complex models and big data sets. For example, 8 to 10+ chains can be run, ensuring convergence quickly by running chains across nodes.

Update I later found it was much easier to run via EasyBuild with OpenMPI using the command runmpi.

Interesting, so this might be an MPI permissions issue then. I will run this by the university IT folks and see if they can help me figure this out. I'll report back.

This seems to most likely have been an issue with mpi on my system, rather than with MrBayes, and this now works for me, although I'm not sure what changed.


text:
How do I interpret the function O(f^-2 log n) in Shekhar et al. (2017)?
The results of Shekhar et al. (2017) "indicate that ASTRAL requires
$\mathcal{O}(f^{-2} \log{n})$

gene trees to reconstruct the species tree correctly with high probability where n is the number of species and f is the length of the shortest branch in the species tree.
How do I interpret this function, specifically the big O notation, when applying it to my data to estimate whether I have enough loci?

Please [edit] your question and provide a link or a full citation to the paper you are quoting. There are 86 papers in PubMed matching (Shekhar[Author]) AND ("2017"[Date - Publication] : "2017"[Date - Publication]) . Also clarify what you mean by applying this function. Why would you apply it?

I've updated your post and added a link to the actual paper you mean. A preprint is available here: https://arxiv.org/abs/1704.06831

Big O notation means "of the order of". Thus $\mathcal{O}(f^{-2} \log{n})$ means that the number of trees required goes with $\log{n}$ for number of species (i.e. if you square the number of genes, you should double the number of trees) and with the inverse square of the shortest branch length (i.e. if you double the shortest branch length, you quarter the number of trees needed).
Big O notation generally allows to to compare two input sizes, but cannot really tell you anything about a single point. I.e. it is hard (impossible?) to tell from this function whether you have sufficient loci.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to count the number of mapped read in 100-bp window from a BAM/SAM file
Although I know how to get total number of mapped read using samtools flagstat (samtools flagstat file_sorted.bam) but I want to count total number of mapped read in a non-overlapping sliding window of fixed size (let's say 500 bp) only with respect to reference (hg19) (assigning each read only once by its start position).
An example of the desired output would be something like:
Chromosome    Region_Start     Region_Stop    Num_Unique_Reads
chr2          1                500            200
chr2          501              1000           202
chr2          1001             1501           230
chr2          .                 .              .
chr2          .                 .              .
chr2          .                 .              .
chr2          8001              8500           20

Can anyone suggest me how to perform that window operation?

one-liner
Here's a gritty one-liner to count the number of reads in a region if you have just one region that you want to investigate. Change the 1 in ($4 >=1) and the 500 in ($4 <=500) to set your window. Change "hg19" to your target sequence.
Note, this one-liner does not double-count reads because of uniq.
samtools view file_sorted.bam | \
    awk '{ if (($3 == "hg19") && ($4 >=1) && ($4 <=500)) { \
  print($1) }}' | sort | uniq | wc -l

python 3
If you need something a little more programmatic to get a known region you should consider trying the pysam package.
import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
region_set = set()
start = 100
stop =  500
counter = 0
for read in samfile.fetch('interesting_contig', 100, 500):
    region_set.add(read.query_name)
    counter += 1
print("{} unique reads in [{}, {}]".format(len(region_set), start, stop))
print("{} alignments in [{}, {}]".format(counter, start, stop))

For my bam file, this program outputs:
1618 unique reads in [100, 500]
2487 alignments in [100, 500]

python 3, non-overlapping windows
It is important to note that the samfile.fetch() function conforms to 1-based indexing, consistent with SAM files, unlike pysam's otherwise default 0-based indexing that matches the BAM-indexing.
import pysam
samfile = pysam.AlignmentFile("file_sorted.bam", "rb")
print(samfile.references)
print(samfile.lengths)
window = 500
print("ref\tstart\tstop\tuniq_reads")
for i in range(len(samfile.references)):
    refname = samfile.references[i]
    seqlen = samfile.lengths[i]
    for j in range(1, seqlen, window):
        stop = j+window-1 if j+window-1 < samfile.lengths[i] else samfile.lengths[i]
        region_set = set()
        for read in samfile.fetch(refname, j, stop):
            region_set.add(read.query_name)
        print("{}\t{}\t{}\t{}".format(refname, j, stop, len(region_set)))

For the small bam file I looked at the output was:
ref         start  stop  uniq_reads
OdonB       0      500   1668
OdonB       501    1000  2647
OdonB       1001   1500  129400
OdonB       1501   2000  598
OdonB       2001   2500  538
OdonB       2501   3000  872
OdonB       3001   3500  561
OdonB       3501   4000  49120
OdonB       4001   4500  48433
OdonB       4501   4806  28121

Thanks for  your reply with proper examples. But, it seems I have missed something in question. By a window, I mean that it is a sliding non-overlapping window of width 500 bp (not just one region). Once again thanks for your answer. Can you extend your answer for sliding non-overlapping windows of given width.

This updated answer seems like it is what you are actually after, right?

Thanks. Can I call number of uniq_reads in a given window as Read Depth of that window (because final results is clearly telling the this many reads are supporting particular region of reference genome.)?

It is pretty similar to "Read Depth", but that makes more sense on a per-position basis. I think it makes sense just to call it what it is -- the number of unique reads mapping to each 500-bp window.

using a streaming process with bioalcidaejdk rather than running many random-access jumps.
The script: the.code
final int WINDOW_SIZE=500;

// map position to count number of reads
final Map<Integer,Long> pos2count = new TreeMap<>();
String currChrom=null;

for(;;)
    {
    // get next read in stream
    final SAMRecord rec = iter.hasNext()?iter.next():null;
    //ignore 'bad' reads
    if(rec!=null ) {
        if(rec.getReadUnmappedFlag()) continue;
        if(rec.getDuplicateReadFlag()) continue;
        if(rec.isSecondaryOrSupplementary()) continue;
        }
    // loop over the windows before the current
    final Iterator<Integer> posIter = pos2count.keySet().iterator();
    while(posIter.hasNext())
        {
        final int start = posIter.next();
        final int end = start + WINDOW_SIZE;
        if(rec!=null && rec.getContig().equals(currChrom) && end >= rec.getStart()) break;
        final long count = pos2count.get(start);
        out.println(currChrom+"\t"+start+"\t"+end+"\t"+count);
        //remove this position from pos2count
        posIter.remove();
        }
    //no more read, exit
    if(rec==null) break;

    //retrieve all the windows for the current read
    for(Integer pos : IntStream.rangeClosed(rec.getStart(), rec.getEnd()).
            mapToObj(POS-> WINDOW_SIZE*(int)(POS/(1.0*WINDOW_SIZE))).
            collect(Collectors.toSet()))
        {
        //increment the count of reads
        final Long count = pos2count.get(pos);
        pos2count.put(pos,1L + (count==null?0L:count));
        }

    currChrom = rec.getContig();
    }

execute:
$java -jar dist/bioalcidaejdk.jar -f the.code  src/test/resources/S1.bam

RF01    0   500 38
RF01    500 1000    61
RF01    1000    1500    69
RF01    1500    2000    80
RF01    2000    2500    73
RF01    2500    3000    56
RF01    3000    3500    20
RF02    0   500 39
RF02    500 1000    76
RF02    1000    1500    73
RF02    1500    2000    72
RF02    2000    2500    59
RF02    2500    3000    13
RF03    0   500 34
RF03    500 1000    74
RF03    1000    1500    85
RF03    1500    2000    67
RF03    2000    2500    51
RF03    2500    3000    7
RF04    0   500 31
RF04    500 1000    68
RF04    1000    1500    90
RF04    1500    2000    75
RF04    2000    2500    29
RF05    0   500 47
RF05    500 1000    100
RF05    1000    1500    53
RF05    1500    2000    5
RF06    0   500 49
RF06    500 1000    89
RF06    1000    1500    29
RF07    0   500 58
RF07    500 1000    65
RF07    1000    1500    7
RF08    0   500 58
RF08    500 1000    69
RF08    1000    1500    3
RF09    0   500 58
RF09    500 1000    69
RF09    1000    1500    2
RF10    0   500 53
RF10    500 1000    39
RF11    0   500 50
RF11    500 1000    32

Thanks. I work in Python and R. So python scripts will better help me. So your approach may help others who work in java.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Snakemake and SLURM -unable to use it with with --cluster
I'm trying to automate some of my previous bash scripts with snakemake. Below is one of the recipe I created:
import os
shell.executable("/bin/bash")

workdir: "/proj/uppstore2018136/publication/jupyterNotebooks"

ENSEMBL_RELEASE='94'
ENSEMBL_GRCm38_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/fasta/mus_musculus/dna'
ENSEMBL_GRCm38_GTF_BASE='ftp://ftp.ensembl.org/pub/release-'+ENSEMBL_RELEASE+'/gtf/mus_musculus'
GTF_FILE='Mus_musculus.GRCm38.'+ENSEMBL_RELEASE+'.gtf'

rule all:
    input:
        "genome/genome.fa",
        "genome/"+GTF_FILE,
        "indexes/bowtie2/reference.1.bt2",
        "indexes/bowtie2/reference.2.bt2",
        "indexes/bowtie2/reference.3.bt2",
        "indexes/bowtie2/reference.4.bt2",
        "indexes/bowtie2/reference.rev1.bt2",
        "indexes/bowtie2/reference.rev2.bt2",
        "indexes/hisat2/reference.1.ht2",
        "indexes/hisat2/reference.2.ht2",
        "indexes/hisat2/reference.3.ht2",
        "indexes/hisat2/reference.4.ht2",
        "indexes/hisat2/reference.5.ht2",
        "indexes/hisat2/reference.6.ht2",
        "indexes/hisat2/reference.7.ht2",
        "indexes/hisat2/reference.8.ht2"

rule getGenome:
    output:
        output1="genome/genome.fa",
        output2="genome/"+GTF_FILE
    shell:
           '''
           curl -o genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz {ENSEMBL_GRCm38_BASE}/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           gunzip genome/Mus_musculus.GRCm38.dna.primary_assembly.fa.gz
           mv genome/Mus_musculus.GRCm38.dna.primary_assembly.fa genome/genome.fa
           curl -o genome/{GTF_FILE}.gz {ENSEMBL_GRCm38_GTF_BASE}/{GTF_FILE}.gz
           gunzip genome/{GTF_FILE}.gz
           '''

rule bowtie2Build:
    input: "genome/genome.fa"
    params:
        basename="indexes/bowtie2/reference"
    output:
        output1="indexes/bowtie2/reference.1.bt2",
        output2="indexes/bowtie2/reference.2.bt2",
        output3="indexes/bowtie2/reference.3.bt2",
        output4="indexes/bowtie2/reference.4.bt2",
        outputrev1="indexes/bowtie2/reference.rev1.bt2",
        outputrev2="indexes/bowtie2/reference.rev2.bt2"
    threads: 20
    shell:
        '''
        module load bowtie2
        bowtie2-build --threads {threads} {input} {params.basename}
        '''

rule hisat2Build:
    input: "genome/genome.fa",
        "genome/"+GTF_FILE
    params:
        basename="indexes/hisat2/reference"
    output:
        output1="indexes/hisat2/reference.1.ht2",
        output2="indexes/hisat2/reference.2.ht2",
        output3="indexes/hisat2/reference.3.ht2",
        output4="indexes/hisat2/reference.4.ht2",
        output5="indexes/hisat2/reference.5.ht2",
        output6="indexes/hisat2/reference.6.ht2",
        output7="indexes/hisat2/reference.7.ht2",
        output8="indexes/hisat2/reference.8.ht2",
    threads: 20
    shell:
        '''
        module load HISAT2
        hisat2_extract_splice_sites.py genome/{GTF_FILE} > genome/genome.ss
        hisat2_extract_exons.py genome/{GTF_FILE} > genome/genome.exon
        hisat2-build -f -p {threads} {input} --ss genome/genome.ss --exon genome/genome.exon {params.basename}
        rm genome/genome.ss
        rm genome/genome.exon
        '''

rule clean:
    shell:
        '''
        rm genome/*
        rm indexes/bowtie2/*
        rm indexes/hisat2/*
        '''

And this is my driver script which I'm submitting as a job itself:
#!/bin/bash

module load snakemake

snakemake --jobs 1 --local-cores 20 --cluster-config /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/cluster.json --latency-wait 120 --max-jobs-per-second 1 -c "sbatch -A {cluster.account} -p {cluster.partition} -n {-cluster.n} -t {cluster.time} -C {cluster.C}" -s indexMake

But I'm getting this error:
Building DAG of jobs...
Using shell: /bin/bash
Provided cluster nodes: 1
Job counts:
    count   jobs
    1   all
    1   bowtie2Build
    1   getGenome
    1   hisat2Build
    4

rule getGenome:
    output: genome/genome.fa, genome/Mus_musculus.GRCm38.94.gtf
    jobid: 2

RuleException in line 30 of /crex/proj/uppstore2018136/publication/jupyterNotebooks/recipes/indexMake:
NameError: The name '-cluster' is unknown in this context. Please make sure that you defined that variable. Also note that braces not used for variable access have to be escaped by repeating them, i.e. {{print $1}

}

I really am not sure what I'm doing wrong here ? Is there a conflict with variables or there is something wrong with the snakemake version I'm using (4.5.0).

(vote to leave open) I know that Snakemake can be used for non-bioinformatics problems. But in practice it is used solely by bioinfo community therefore I think debugging snakefiles is appropriate for this Q&A site.

(vote to close) The official Snakemake docs point folks to Stack Overflow. I suggest we use StackExchange's migrate mechanism to send generic Snakemake questions to SO to keep it all in one place.

I'm voting to close this question as off-topic because it belongs on StackOverflow (although I lack the reputation to initiate the migration).

@DanielStandage it does seem to be aimed primarily at bioinformaticians, and is indeed a tool commonly used by the community, so it seems reasonable that we should support it.

In your submission script you have -n {-cluster.n}, hence the error from snakemake I guess. Are you sure it is not a typo and it should be -n {cluster.n} instead?

Yes, that was indeed the error. I'm not sure how it crept in - working on Sundays is not that productive -


text:
Rdkit.Chem.Draw.savefig() creates empty figures
I can draw a molecule in a jupyter notebook with RDKit. Though, when I use the savefig() method only a white canvas is saved. No molecule is saved.
%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')
fig.savefig('mol.jpeg')

Any idea what goes wrong here?

For some reason you have to set bbox_inches to tight. Without that rdkit does not fit the molecule into the canvas and might plot an empty figure.
fig.savefig('mol.jpeg', bbox_inches='tight')


text:
How to remove frame on RDKit figures?
RDKit can plot molecules, thought the structures are surrounded by an ugly frame with ticks. How can I turn that off and plot the molecules without that frame?
Note, the %pylab inline portion of this code is from the Jupyter environment.
%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')

It looks like rdkit and rdkit's Chem just use the matplotlib. So, if you import matplotlib you can edit your figure like any other matplotlib figure.
I produced this from the command line, so it may work differently in the Jupyter environment.
It looks like Chem doesn't properly scale the drawing based on how large the chemical is, so I added plt.xlim and plt.ylim to make the molecule larger in the frame.
from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
plt.title('Test')
plt.axis('off')
coord_min = 0.35
coord_max = 0.65
plt.xlim((coord_min, coord_max))
plt.ylim((coord_min, coord_max))
fig.savefig('mol.jpeg', bbox_inches='tight')


text:
R package grur: missing_visualization() fails during redundancy analysis of .vcf
I am trying to use the R package grur to process GBS SNP data, which is installed from github as instructed.  For starters, I am running the test code in this vignette.  As suggested I am running it from console in RStudio (v. 1.1.383), and with a recent R installation (below).  I do, and get:
 > ibm <- grur::missing_visualization(data = "example_vcf2dadi_ferchaud_2015.vcf", strata = "strata.stickleback.tsv")
#######################################################################
#################### grur::missing_visualization ######################
#######################################################################
Folder created:
missing_visualization_20181209@1523

Importing data

Reading VCF...
Generated a filters parameters file: filters_parameters_20181209@1523.tsv

Number of SNPs: 31802
Number of samples: 177

conversion timing: 1 sec

VCF: biallelic SNPs
Cleaning VCF sample names

Synchronizing sample IDs in VCF and strata...
Reads assembly: reference-assisted
Filters parameters file: updated

Number of chromosome/contig/scaffold: 196
Number of locus: 17095
Number of markers: 31802
Number of individuals: 177

Working time: 1 sec

Deprecated function, update your code to use: filter_monomorphic

Scanning for monomorphic markers...
Number of markers before/blacklisted/after: 31802/0/31802

Tidy genomic data:
    Number of markers: 31802
    Number of chromosome/contig/scaffold: 196
    Number of individuals: 177
    Number of populations: 8

Informations:
Number of populations: 8
Number of individuals: 177
Number of ind/pop:
HAD = 21
HAL = 19
KIB = 17
KRO = 20
MOS = 20
MAR = 20
NOR = 20
ODD = 40

Number of duplicate id: 0
Number of chrom/scaffolds: 196
Number of locus: 17095
Number of SNPs: 31802

Proportion of missing genotypes (overall): 0.022169

Identity-by-missingness (IBM) analysis using
    Principal Coordinate Analysis (PCoA)...
Generating Identity by missingness plot

Redundancy analysis...

Error in UseMethod("ungroup") :
  no applicable method for 'ungroup' applied to an object of class "character"
In addition: Warning message:
Trying to compute distinct() for variables not found in the data:
- `strata.select`
This is an error, but only a warning is raised for compatibility reasons.
The operation will return the input unchanged.

The final object is not generated, however some plots of the PCoA analysis are. I get this result on two Macs, one with Sierra, the other with High Sierra, but both with the same RStudio and R versions. All R packages are up to date. I get the same error when I run my own .vcf
Does anyone else get these errors on this test dataset from the repo? Can anyone suggest how to deal with this?
> version
platform       x86_64-apple-darwin15.6.0
arch           x86_64
os             darwin15.6.0
system         x86_64, darwin15.6.0
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

This seems to be an internal error (using dplyr::ungroup without doing dplyr::group) . If you provide a minimal reproducible example, the maintainer of the package might be able to solve it (The maintainer should be able to reproduce the error, so provide enough data to make it possible).


text:
Importing GFF file with Biopython
Is there a way to import a GFF file for an organism with Biopython in the same way you can for a Genbank file?
For example,
from Bio import Entrez as ez
ez.email = '...'
handle = ez.efetch(db='gene', id='AE015451.2', rettype='genbank', retmode='xml')
h = handle.read()
print(h)

This will import a Gebank xml file for Pseudomonas putida.  Is there a way to do this for a GFF file?

No, there is currently no GFF support in biopython.
However, you can read in GFF files into python using this package, gffutils. There are also a few other packages to read/write GFF files, like gff3.

I would use BCBio for gff handling as it is written to directly interface with BioPython’s object model.
The only downside is that I believe it is no longer actively supported. It is however the package that the BioPython docs use generally. There are plans to properly incorporate GFF/GTF parsers in to BP in the not too distant future according to that link though.

Most GFF parsers handle the work of reading annotation data into objects for convenient data access, but most do not handle the important task of resolving relationships between features. The tag Python library does both, grouping related features together for inspection, traversal, and feature-by-feature processing.
CAVEAT: The tag library only supports GFF3.
DISCLAIMER: I'm the developer of the tag library, so this is a shameless plug. :-)


text:
How BLOSUM Matrix is constructed and calculated
I would like to ask how BLOSUM matrix is constructed and calculated ? I read the wikipedia and I still do not understand it. I do not understand the mathematical calculations as I have low knowledge about logarithms.
Maybe could someone give a simple example and then a real-scenario example of how the BLOSUM matrix can be used and calculated, maybe in relation to the Blast algorithm where it uses the BLOSUM scoring matrix to determine high scoring words for each word in the query sequence ?
Also, why are they different types of BLOSUM matrices like BLOSUM62 and so on ?

Sean Eddy wrote an excellent explanation of the BLOSUM Matrix in his Nature Biotechnology paper Where did the BLOSUM62 alignment score matrix come from? He also gives real examples of how one would calculate the score for different amino acid pairs and where the numbers come from.  Further references given at the end of the paper might also be useful reading, specifically: Amino acid substitution matrices from an information theoretic perspective, Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes, and Amino acid substitution matrices from protein blocks.

Could you please add a summary of what these links say to your answer? Answers need to be able to stand on their own even if the inks are removed.

Hi @Matt Bashton, I can't access the reading materials in the links as It needs subscription to read it.

BLOSUM stands for BLOcks SUbstitution Matrices (Henikoff & Henikoff, 1992), and were created by observing substitution frequencies in local ungapped multiple sequence alignments. The score reflects the chance (log-odds) one amino acid is substituted for another in a set of protein multiple sequence alignments. BLOSUM scoring matrices are normally followed by a number eg BLOSUM62. The number is % sequence identity between the sequences in the multiple sequence alignment (MSA) used to create the score matrix. eg BLOSUM62 matrices were created from multiple sequence alignments with blocks that shared 62% identity. Note because these are ungapped alignments, the multiple sequence alignments were generated from conserved protein families.  With the older PAM matrices, the number is opposite, a higher number is greater distance. PAM matrices reflect evolutionarily distant in global (aka gapped) alignments of the full-length sequences that include both conserved and diverged regions.  See Wikipedia https://en.wikipedia.org/wiki/BLOSUM. There are several newer substitution matrices.  The R package SubVis allows you to compare/visualize different substitution matrices.

Hi @Aedin, thanks for your input. Is it possible to provide an example to run through the formula itself and also showing how the matrix is formed as shown in the Blosum matrix in the wikipedia ?

Amino acid matrices are contructed by maximum likelihood and easily performed using PAML. You present a known tree and an amino acid alignment and ask what amino acid matrix would maximise the likelihood. Hence all amino acid matrices using this old school approach are dependent on the alignment composition. JTT was generally regarded as obtaining the highest likelihood. However, I reitereate that things have moved on a lot since this 2-D reversible amino acid substitution matrix was devised.


text:
Derive a GTF containing protein coding genes from a GTF file with Exons and CDS
Why I need a compatible file
I’m trying to run
velocyto with the R package to analyse RNA velocity (cell trajectories) with single cell RNASeq data. I have performed single cell analysis from 10x Genomics data using cellranger.
I have successfully aligned the reads to get loom files and imported these into R. I can get the velocity from these files by following the vignettes. However, I cannot reproduce the RNA velocity analysis based on “gene structure”.
I’m working with a different organism to the example (not human or mouse) so the annotation data provided in the example does not work. I have a GTF file for the latest annotation of this organism. However, it only contains “exon” and “CDS” as features. This appears to be the source of the problem. The “find.ip.sites” function requires a GTF with “features” = “gene” and one of the “attributes” to be “protein_coding”. These requirements are hard-coded into the velocyto.R function.
I have the following GTF files from the AtRTD2 dataset. The chromosome labels match the loom files in R.
Chr1    TAIR10  exon    3631    3913    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    3996    4276    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4486    4605    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    4706    5095    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5174    5326    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  exon    5439    5899    .   +   .   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3760    3913    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 3996    4276    .   +   2   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4486    4605    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 4706    5095    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5174    5326    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    TAIR10  CDS 5439    5630    .   +   0   transcript_id "AT1G01010.1"; gene_id "AT1G01010"; gene_name "AT1G01010";
Chr1    Araport11   exon    6788    7069    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7157    7450    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7564    7649    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7762    7835    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    7942    7987    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8236    8325    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8417    8464    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   exon    8571    8737    .   -   .   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7315    7450    .   -   1   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7564    7649    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7762    7835    .   -   2   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 7942    7987    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";
Chr1    Araport11   CDS 8236    8325    .   -   0   transcript_id "AT1G01020_P2"; gene_id "AT1G01020"; gene_name "AT1G01020"; Note "ARV1";

The specifications for a GTF or GFF2 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqname - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
feature - feature type name, e.g. Gene, Variation, Similarity
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
frame - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attribute - A semicolon-separated list of tag-value pairs, providing additional information about each feature.

What I’m looking for
Is there a way to derive a compatible GTF file containing protein coding genes from the exons and CDS? I want to produce a GTF which contains genes in the features and protein_coding in the attributes. Is it possible to do this with existing tools or scripts?

What I’ve tried so far
I can modify the source code of the “find.ip.sites” function to run on my GTF file with these features missing. However, this requires running internal functions from the package written in Rcpp and means my workflow override future updates to the package. Running the rest of the vignette returns errors as no introns long enough have been identified (setting the thresholds lower is also incompatible with the general linear models called). Therefore I think it is better to generate a compatible GTF or GFF3 file rather than alerting the source code of the functions.
While intended for GTF files, the package functions can import GFF3 files, despite being described for GTF in the documentation. I’ve tried generating a GFF3 file using gffread from Cufflinks and replacing “feature” = “mRNA” with “gene” and adding “protein_coding”. This also returns errors when running the velocity algorithm. It does not work for either the GTF file used as an input for cellranger or the file generated by it. Is there a way to annotate protein coding genes and intron/exon boundaries based on a GTF file containing only exon and CDS annotations?
With cufflinks 2.2.1, a GFF3 was generated
gffread - E AtRTD2_19April2016.gtf -o- > AtRTD2_19April2016.gff
sed -i '/mRNA/s/gene_name=AT/gene_type=protein_coding;gene_name=AT/g' AtRTD2_19April2016.gff'
sed -i 's/mRNA/gene/g' AtRTD2_19April2016.gff

This is the GFF3 files that I've tried:
# gffread - E AtRTD2_19April2016.gtf -o-
##gff-version3
Chr1    TAIR10  gene    3631    5899    .   +   .   ID=AT1G01010.1;geneID=AT1G01010;gene_type=protein_coding;gene_name=AT1G01010
Chr1    TAIR10  exon    3631    3913    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    3996    4276    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4486    4605    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    4706    5095    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5174    5326    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  exon    5439    5899    .   +   .   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3760    3913    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 3996    4276    .   +   2   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4486    4605    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 4706    5095    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5174    5326    .   +   0   Parent=AT1G01010.1
Chr1    TAIR10  CDS 5439    5630    .   +   0   Parent=AT1G01010.1
Chr1    Araport11   gene    6788    8737    .   +   .   ID=AT1G01020_P2;geneID=AT1G01020;gene_type=protein_coding;gene_name=AT1G01020
Chr1    Araport11   exon    6788    7069    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7157    7450    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7564    7649    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7762    7835    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    7942    7987    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8236    8325    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8417    8464    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   exon    8571    8737    .   -   .   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7315    7450    .   -   1   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7564    7649    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7762    7835    .   -   2   Parent=AT1G01020_P2
Chr1    Araport11   CDS 7942    7987    .   -   0   Parent=AT1G01020_P2
Chr1    Araport11   CDS 8236    8325    .   -   0   Parent=AT1G01020_P2

The specifications for a GTF or GFF2 are:
The specifications for a GF3 are:

Fields
Fields must be tab-separated. Also, all but the final field in each feature line must contain a value; "empty" columns should be denoted with a '.'

seqid - name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seq ID must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
source - name of the program that generated this feature, or the data source (database or project name)
type - type of feature. Must be a term or accession from the SOFA sequence ontology
start - Start position of the feature, with sequence numbering starting at 1.
end - End position of the feature, with sequence numbering starting at 1.
score - A floating point value.
strand - defined as + (forward) or - (reverse).
phase - One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
attributes - A semicolon-separated list of tag-value pairs, providing additional information about each feature. Some of these tags are predefined, e.g. ID, Name, Alias, Parent - see the GFF documentation for more details.

It would probably help of you could [edit] your question and show us a few lines of each of the two GTF files so we can see what you mean more clearly.

You probably can use AWK for this, is the "gene" feature similar as "CDS" or is it from "first exon start" to "last exon end".

I think the question and the "What I'm looking for" are a little confusing. Specifically, it sounds like you want to generate a GTF in "What I'm looking for" but you want to get the amino acid sequence of proteins from a GTF in the title. Is there any way you can clarify or make the title and "What I'm looking for" consistent?

We have a complete GTF for Arabidopsis in Ensembl Plants.

Thanks @Emily_Ensembl, I was able to get velocyto to compute "global" RNA velocity based on this reference sequence. In case anyone else runs into this issue, the loom file needs to be re-computed with same GTF file with the shell utility to be compatible with the R commands.

In your case, I would definitely suggest following @Emily_Ensembl's advice and using the Arabidopsis GTF from Ensembl. But for future reference, if an Ensembl GTF wasn't available, you could build something like this using the gtf class from cgat
The cgat module and dependancies can be installed by following the instructions here. Specifically, you can download and run their installation script. This will set up a dedicated conda environment to run the versions of dependancies needed. Further documentation can be found here.
# download installation script:
curl -O https://raw.githubusercontent.com/cgat-developers/cgat-apps/master/install.sh

# install the development version (recommended, no production version yet):
bash install.sh --devel

This requires python3 and anaconda to install. You may need to run this script in the conda environment:
# show available environments
conda info --envs
# activate by path
source activate /home/user/local/bin/conda-install/envs/cgat-a

Then run the follow as a python3 script. Copy the contents in a file such as convert_gtf.py. Then run it in the terminal on your gtf file:
python convert_gtf.py genes.gtf > genes_new.gtf

#python 3.6.4

import sys

from cgat import GTF
from cgatcore import iotools

for gene in GTF.flat_gene_iterator(
     GTF.iterator(iotools.open_file(sys.argv[1])), strict = False):
    gene_start = min (e.start for e in gene)
    gene_end = max(e.end for e in gene)

    if any(e.feature == "CDS" for e in gene):
        gene_type = "protein_coding"
    else:
        gene_type = "non_coding"

    for exon in gene:
        exon_line.setAttribute("gene_biotype", gene_type)

    gene_line = GTF.Entry().fromGTF(gene[0])
    gene_line.feature = "gene"
    gene_line.start = gene_start
    gene_line.end = gene_end

    # Its not clear what the transcript_id for a gene line is.
    # Technically speaking, Ensembl GTFs do not follow the GTF
    # standard as GTF entries must include transcript_ids, but
    # the ensembl ones don't.

    if hasattr(gene_line, "gene_name"):
        gene_line.attributes["gene_name"] = gene[0].gene_name
    else:
        gene_line.attributes["gene_name"] = gene[0].gene_id

    gene_line.attributes["gene_biotype"] = gene_type

    print (str(gene_line)+"\n")
    for exon in gene:
        print(str(exon_line)+"\n")

That’s great to know. Since this pertains to using tools in R, it would help to have instructions to install the library and run for those unfamiliar with Python.

I've tested this script and made minor changes. This is what worked for me (to covert the GTF: I'm not sure if it's compatible with velocyto yet).

Thanks for the fixes. The things have changed recently and I lose track of what version things are on. I believe the CGAT team is working on simplifying the installation. A couple of points: I'm not quite sure what you are trying to achieve with the if hasattr(gene_line, "gene_name) bit: presumably if it already has a gene_name attribute, you don't need to set one? Also note that your resulting lines arn't going to have a transcript_id. While this is also true of Ensembl GTFs, it does mean these GTFs don't mean the  GTF specification.

I had issues with missing gene names and often TAIR IDs are used instead of gene symbols. So of this might be specific to Arabidopsis. Most important is the fixes to the loop for exons as these lines weren’t defined.

Its not that I don't understand that you might want to add a gene_name to the attributes, its just that currently you test for the presence of a gene_name and then set it to its already currently defined value if it already has one. This seems redundant.


text:
What should I use as a reference file for clumping?
I'm new to using PLINK and I have 2 basic questions about clumping. I'm computing polygenic risk scores for psychiatric diseases and will be using the following command to clump:
plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb

My first question: What type of file should I be using for the reference? Is it a map file?
My second question: How do I choose the most appropriate reference file?


text:
Snakemake: Cannot find first rule?
Back again with another snakemake query. This time I decided to port my read cleaning and alignment pipeline to snakemake. Repeating the steps from previous question and trying not to make any typo here, I'm now receiving this "error" while trying a dry run using snakemake snakemake -nr -s alignMake:
Building DAG of jobs...
MissingRuleException:
No rule to produce snakemake (if you use input functions make sure that they don't raise unexpected exceptions).

Below is my snakefile:
import os
shell.executable("/bin/bash")

workdir: '/proj/uppstore2018136/publication/jupyterNotebooks'

pathToPEReads = 'raw_data/reads/BL6/totalRNA'
pathToSMReads = 'raw_data/reads/BL6/smallRNA'

pathToHISAT2index = 'indexes/hisat2/reference'
pathToBOWTIE2index = 'indexes/bowtie2/reference'

pathToTemp = 'temp/'
pathToBAM = 'processed_data/BAM'

readsPE, = glob_wildcards('raw_data/reads/BL6/totalRNA/{sample}_1.fastq.gz')
readsSE, = glob_wildcards('raw_data/reads/BL6/smallRNA/{sample}.fq.gz')

samToolsThreads = 6

rule all:
    input:
        expand('processed_data/BAM/{sample}.bam', sample=readsPE),
        expand('processed_data/BAM/{sample}.bam', sample=readsSE)

rule trimgalorePE:
    input:
        fastq1 = pathToPEReads + '{sample}_1.fastq.gz',
        fastq2 = pathToPEReads + '{sample}_2.fastq.gz'
    output:
        output1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        output2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --paired --length 80 --illumina --stringency 2 --clip_R1 11 --clip_R2 11 -o {pathToTemp} --paired {input.fastq1} {input.fastq2}
        '''

rule trimgaloreSEandDust:
    input:
        fastq = pathToSMReads + '{sample}.fq.gz',
    output:
        output = pathToTemp + '{sample}_trimmed.fq.gz',
    shell:
        '''
        module load bioinfo-tools
        module load cutadapt
        module load TrimGalore
        trim_galore --quality 20 --length 16 --max_length 32 --small_rna --fastqc --stringency 2 -o {pathToTemp} {input.fastq}
        '''

rule align_hisat2:
    input:
        fastq1 = pathToTemp + '{sample}_1_val_1.fq.gz',
        fastq2 = pathToTemp + '{sample}_2_val_2.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load HISAT2
        module load samtools
        hisat2 -q -p {threads} --no-discordant --no-mixed --rna-strandness RF --dta -x {pathToHISAT2index} -1 {input.fastq1} -2 {input.fastq2} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule align_bowtie2:
    input:
        fastq = pathToTemp + '{sample}_trimmed.fq.gz'
    output:
        output = 'processed_data/BAM/{sample}.bam'
    threads: 10
    shell:
        '''
        module load bioinfo-tools
        module load bowtie2
        module load samtools
        bowtie2 --no-unal --end-to-end -D 20 -R 3 -N 1 -L 20 -i S,1,0.50 -p {threads} -x {pathToBOWTIE2index} -U {fastq} | samtools sort -@ {samToolsThreads} -o {output}
        '''

rule clean:
    shell:
        '''
        rm temp/*
        '''

I believe it is due snakemake being unable to find the first rule ? Maybe there is something wrong with parsing ? I tested the glob_wildcards and expand function independently and they appear to work properly.
Edit: Here is the directory tree:
.
├── citations
├── doc
├── genome
│   ├── genome.fa
│   └── Mus_musculus.GRCm38.94.gtf
├── indexes
│   ├── bowtie2
│   │   ├── reference.1.bt2
│   │   ├── reference.2.bt2
│   │   ├── reference.3.bt2
│   │   ├── reference.4.bt2
│   │   ├── reference.rev.1.bt2
│   │   └── reference.rev.2.bt2
│   └── hisat2
│       ├── reference.1.ht2
│       ├── reference.2.ht2
│       ├── reference.3.ht2
│       ├── reference.4.ht2
│       ├── reference.5.ht2
│       ├── reference.6.ht2
│       ├── reference.7.ht2
│       └── reference.8.ht2
├── LICENSE
├── notebooks
├── plots
├── processed_data
│   └── BAM
├── raw_data
│   └── reads
│       └── BL6
│           ├── smallRNA
│           │   ├── smBL6_1.fq.gz
│           │   ├── smBL6_2.fq.gz
│           │   └── smBL6_3.fq.gz
│           └── totalRNA
│               ├── EKIxxxx012_1.fastq.gz
│               ├── EKIxxxx012_2.fastq.gz
│               ├── EKIxxxx013_1.fastq.gz
│               ├── EKIxxxx013_2.fastq.gz
│               ├── EKIxxxx014_1.fastq.gz
│               ├── EKIxxxx014_2.fastq.gz
│               ├── ESSxxxx077_1.fastq.gz
│               └── ESSxxxx077_2.fastq.gz
├── README.md
├── recipes
│   ├── alignMake
│   ├── cluster.json
│   ├── indexMake
│   ├── indexMakeScript.out
│   └── indexMakeScript.sh
├── results
├── scripts
├── slurmLogs
│   ├── slurm-6575541.out
│   ├── slurm-6575543.out
│   ├── slurm-6575815.out
│   └── slurm-6580867.out
├── source
└── temp

Loosely related, I think your rule clean is never executed unless your run snakemake as snakemake [opts] clean. If you want to clean up stuff on completion, maybe look at the onsuccess and similar handlers. Also consider the temp flag...

@dariober - I agree. I didn't knew about onsuccess. I will include such handlers in future scripts. But the reason why I didn't want to automatically use clean was to manually check if things happened properly in background.

you are starting snakemake the wrong way.
snakemake snakemake

means "Running the pipeline in Snakefile until you reach a rule named snakemake or a file name snakemake is created."
If your snakefile is called snakemake and this is the pipeline you want to start, the correct syntax would be:
snakemake -s snakemake

Thanks - another one of the stupid mistake - since dry runs are also automated; I didn't notice that the driver was using using snakemake snakemake - should have noticed it while posting the question. Although I discovered it earlier, I guess this question should still stay online.


text:
Genomics: Downloading genomic data from the NCBI Genbank database
I'm taking a genomics course at school, and we have a final project due in a few weeks.
Now, I can write bash scripts to sort and analyze the data, etc. However, I'm not quite sure how to actually get the data from the database. I was hoping to get some help there.
I'd like to get the assemblies of genomes for the species in the Borrelia genus. Is there any way that I can do this using commands like "wget", for example. I'm not even sure what URL I'd have to give the computer to go fetch these. Then once I give it the correct URL, how should I tell it what file in that directory to download?
I'd greatly appreciate some help on this issue.

Have you tried Google on something like Borrelia genome?  I've never done it before but it took me less than 10 minutes to find the genome for a borrelia species from that point.

@mgkrebbs - I'm sorry, but I forgot to mention in my post, that I need the genomes for many species of Borrelia. I wanted to know that, instead of doing it manually, is there some way I could get the computer to download the assemblies for all the Borrelia species in Genbank?

Entrez is what you want for this, most likely. There are R/python packages that can also access it pretty easily (biomaRt/biopython, respectively).

Do you need the assembly information and metadata? or do you need the genomes sequences themselves?

@mgkrebbs Getting data off of NCBI is infuriatingly difficult and non-obvious, if you don’t know exactly what you’re looking for. Finding the data is one thing (already not quite as easy as you make it sound). And then …?!

For downloading data from commandline use the E-utilities of NCBI.

Can you edit your answer to provide an example query to download all of the Borrelia genomes from NCBI using the tool that you linked?

@conchoecia I answered this in [biology.se] and therefore did not go into details. I'll provide an example as soon as I can.

@WYSIWYG could you provide a link to your answer in Biology?


text:
How to plot p-values in a circular barplot?
Im trying to plot reactome pathway in for of circular bar plot using ggplot2
my data
    dput(df2)
structure(list(ID = 1:10, Pathway_names = structure(c(5L, 1L,
7L, 6L, 2L, 3L, 9L, 8L, 10L, 4L), .Label = c("Antigen Presentation: Folding assembly and peptide loading of class I MHC",
"Antigen processing-Cross presentation", "Class I MHC mediated antigen processing & presentation",
"Cytokine Signaling in Immune system", "Endosomal/Vacuolar pathway",
"ER-Phagosome pathway", "Immunoregulatory interactions between a Lymphoid and a non-Lymphoid cell",
"Interferon alpha/beta signaling", "Interferon gamma signaling",
"Interferon Signaling"), class = "factor"), Entities_found = structure(c(1L,
2L, 3L, 4L, 4L, 5L, 6L, 7L, 8L, 9L), .Label = c("51", "54", "55",
"57", "75", "77", "96", "117", "176"), class = "factor"), Entities_pValue = c(1.11e-16,
1.11e-16, 0.003908431, 2.54e-12, 2.77e-10, 0.0053068, 1.64e-13,
1.11e-16, 1.11e-16, 3.68e-06)), row.names = c(NA, -10L), class = "data.frame")

my code
ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found)) + ######################################
  geom_bar(aes(fill = Entities_pValue), stat = "identity", position = "stack") +
  #scale_fill_manual(values = c(blank = "white", dark = "black")) +
  coord_polar() +
  theme_minimal() +
  guides(fill = FALSE)

So far im not sure how to include both p value and entities in the plot so that the pathway that are coming up can be labeled accordingly
New figure

I don't understand the issue. It looks fine to me.

What do you mean by entities: the third column of your df2? You already have the height for the entities and the color for the p-value. How is this plot not what you expected?

yes i tried different ways but this is what i end up with  ,but can you suggest me a way to color it so that it can distinguish p value it color scale_manual ?

@krushnachChandra Sorry I don't understand your question from the comments. You want to distinguish the p-value with the colors? Or do you want to set a manual scale of colors? If the question is the former, you can try to apply log10  to the variable to better distinguish the values. If the later case, I don't understand what is your problem with scale_manual

Can I suggest using acronyms for the entity names in the figure? Like ERPP for ER-Phagosome pathway? This way the text can be placed next to its barplot component and be readable.

okay let me try that

The adjustcolor function from grDevices can be used to generate colours with transparency. What you want is colours based on the pathway and transparency to form a colour gradient for the p-value.
adjustcolor( "red", alpha.f = 0.2)
For the ggplot example, the aes function has an alpha parameter.
ggplot(df2, aes(x = factor(Pathway_names), y = Entities_found), alpha = Entities_pValue) + ######################################
geom_bar(aes(fill = Entities_found, alpha = Entities_pValue), stat = "identity", position = "stack") +
coord_polar() +
theme_minimal() +
 guides(fill = FALSE)

it looks fine now may be i can tweak it a bit more i will see but i see my Entities_found are still 150 but highest one is 176 why am not getting till 176?

i have a doubt the p value gradient still looks confusing ...can you help me regarding that

Try plotting them on a log-scale?

no i haven;t given that thought may be that would help


text:
How to find all variable-length seqs with an exact 5' and 3' match in a FASTA file
Context
I am interested in finding all of the promotors specific to a particular sigma factor. I have identified the -35 and -10 sites from the literature, bold denotes -10, -35, binding sites:
  -35 site               -10 site
5' TTTACAtatttatttcagacaacGTCTTT 3'
   ^^^^^^                 ^^^^^^

The CAPITALIZED and ^emphasized^ 5' and 3' sequence components are the conserved sequence that must have an exact match. The lower-case bases are the variable 10-25bp region for which the sequence does not matter.
The sigma factor is from C. difficile and I am particularly interested in searching C. acetobutylicum, however, it would also be interesting to see where else this sequence crops up. Obviously such a short sequence will arise often by random chance, so it is important to define the spacer region, which might be anywhere between 15-25 nucleotides and may be any sequence.
In particular, I would like to find these binding regions from this sequence and this one.
Desired output
A text file that contains the locations and sequences of all of the variable-width sequences as defined above. Something like:
seq   start  stop  gap  dir  seq
seq1  1      22    10   F    TTTACAagtcagatctGTCTTT
seq2  2      23    11   F    TTTACAagtcagactctGTCTTT
seq2  122    101   10   R    AAAGACatgctagccaTGTAAA

Thank you.

Is this the same as your other question "BLAST Promotor -35 to -10"? If so, please consider to close/delete the other one. Thanks.

python 3 solution
This python program will get you what I think you want. It is not elegant given that it (A) stores the whole sequence in memory and (B) spends most of its time slicing strings rather than working with a buffer. This program is not recommended for FASTA files for which the largest sequence will not fit well in RAM.
All that you need to do to modify it for other genomes or other binding sites is to change the three lines that begin with forseq, revseq, and filename.
from Bio import SeqIO
from Bio.Seq import Seq

forseq = Seq("tttaca".upper())
revseq = Seq("gtcttt".upper())
filename = "NC_003030.fasta"

forcom = forseq.reverse_complement()
revcom = revseq.reverse_complement()
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if record.seq[i:i+len(v1)] == v1:
                for j in range(10,25):
                    if record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] == v2:
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here are the first few lines of example output from the first fasta file that you linked.
seq          start    stop     gap  dir  seq
NC_003030.1  97512    97543    20   F    TTTACAAAAGAAAAGATATGCATAATGTCTTT
NC_003030.1  283173   283205   21   F    TTTACACTTATGTTTATATAGTTCAACGTCTTT
NC_003030.1  303104   303125   10   R    AAAGACGCAAGCATAATGTAAA
NC_003030.1  1127691  1127725  23   F    TTTACAAGGTGTTACTGATAGAACAAAGGGTCTTT
NC_003030.1  1301958  1301988  19   F    TTTACATAAAAAGATGATAGTTTTTGTCTTT

python 3 with redundancy
If you need sequence redundancy, this version works. I used the redundancy codes available from www.reverse-complement.com You may use the redundant characters in [ R, Y, S, W, K, M, B, V, D, H, N ] or lowercase variants. I'm not sure how to handle the '-' character at the moment but am open to suggestions. Code and results below.
Again, the algorithm isn't great given all the for loops so other answers are appreciated!
from Bio.Seq import Seq

forseq = "tttaca" # Just use python strings here. No Seq() as above
revseq = "ntnntn"
filename = "NC_003030.fasta"

revcomp = {'A':'T', 'C':'G', 'G':'C', 'T':'A',
           'R':'Y', 'Y':'R', 'S':'S', 'W':'W',
           'K':'M', 'M':'K', 'B':'V', 'V':'B',
           'D':'H', 'H':'D', 'N':'N'}
redund = {'A':['A'], 'C':['C'], 'G':['G'], 'T':['T'],
          'R':['A','G'], 'Y':['C','T'],
          'S':['C','G'], 'W':['A','T'],
          'K':['G','T'], 'M':['A','C'],
          'B':['C','G','T'], 'V':['A','C','G'],
          'D':['A','G','T'], 'H':['A','C','T'],
          'N':['A','C','G','T']}

def revcom(seq):
    return "".join(revcomp.get(base, base) for base in reversed(seq))

def seqmatch(query, reference):
    """Checks if the query sequence (perhaps a promoter with redundant bases)
    matches the reference. Both query and reference must be the same length.

    Returns true if there is a match (redundancy allowed), returns false
    otherwise.
    """
    assert len(query) == len(reference)
    for i in range(len(query)):
        if reference[i] in redund[query[i]]:
            pass
        else:
            return False
    return True

forseq = forseq.upper()
revseq = revseq.upper()
forcom = revcom(forseq)
revcom = revcom(revseq)
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if seqmatch( v1, record.seq[i:i+len(v1)] ):
                for j in range(10,25):
                    if seqmatch( v2, record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] ):
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here is some sample output:
seq          start  stop  gap  dir  seq
NC_003030.1  733    767   23   R    GATGAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  736    767   20   R    GAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  739    767   17   R    GATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  742    767   14   R    CAAGAAACCGATACAAACAATGTAAA
NC_003030.1  743    767   13   R    AAGAAACCGATACAAACAATGTAAA
NC_003030.1  2665   2697  21   R    TATTAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  2668   2697  18   R    TAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  3798   3822  13   R    CACCAGAAGACTTAAAAATTGTAAA
NC_003030.1  3801   3822  10   R    CAGAAGACTTAAAAATTGTAAA

Thank you so much for this. I really appreciate it. I am however having some issues running the script:

I have installed biopython. Do I need to download 'NC_003030.fasta' and specify the directory?

Yes, you need to change 'NC_003030.fasta' to the path of your file, such as filename = "/I/saved/my/file/here/NC_003030.fasta"

Hello, someone would like to use this code too so I thought I had better check with you for permission to share - even though this is an open source forum it seems like the proper thing to do. Also, I will almost certainly use this in my thesis and will need to reference you. Can you provide a name/affiliation?

Thanks! That is awfully considerate of you to ask. If you would like to do that, while totally not necessary, I found a few relevant SE questions: academia.stackexchange.com/questions/40325/ and meta.stackexchange.com/questions/191943 . A quick search of my username on some other sites popular with scientists/programmers and you may find my contact info. I'm also working on my thesis and usually read 100+ SE questions per day - this is a good thing for me to consider too!

No problem. I guess since you are a dry lab scientist it's very common for you to share scripts with your colleagues without thinking much of it, it would be unrealistic to trace and record the origin of every line of code perhaps. I am a wet lab scientist and I'm not really familiar with the customs amongst working coders. I will be more or less pasting the code into my M&M so it really wouldn't be right to not reference it in my case.

I do have one more question. Only two bases appear to be conserved in the -10 sequence. How might I specify a variable base. I have tried using an 'x' an 'n' or a '-' eg 'ntnntn' none of these variations produce a result, but I do not receive an error message either.
'

Holy nested for loop!

Yep. It's a beast!

@Ryan_J_Hope many are both wetlab and bioinformaticians. You're doing a great thing by taking charge of both research fronts! If you would like an expanded answer allowing the gap character '-' maybe update the question? A version allowing redundancy but no gaps is above in the latest edits.


text:
Albacore basecalling running but outputs 0 reads
I am trying to basecall data produced by the MinION using the SQK-LSK109 kit and FLO-MIN106 flowcell via the command-line. My version of albacore is the latest (v2.3.4). I tried running using the following command-line (the name of the sequencing run and fast5s are different and are default from machine output):
python3.6 /home/dkennetz/.local/bin/read_fast5_basecaller.py -i ./seqrun/fast5/ -t 1 -s ./good_reads/ -f FLO-MIN106 -k SQK-LSK109 -r -o fastq

-t 1 is just for this test. It is running and finishing almost instantly with 0 errors but outputting no reads. I have previously run this same command-line using the SQK-LSK108 kit and had no problems.
I then went on to validate the fast5s using the h5_validate tool:
python3.6 /home/dkennetz/.local/bin/h5_validate /home/dkennetz/.local/lib/python3.6/site-packages/h5_validator/schemas/single_read_fast5.yaml seqrun/fast5/0/seqrun.fast5 -v

The output from this is:
    Matcher <'UniqueGlobalKey', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

Error at /:

Matcher <'PreviousReadInfo', type=group match_type=exact count=
{'minimum_count': 1, 'maximum_count': 1}> was not satisfied after matching <HDF5 file "seqrun.fast5" (mode r)>

I read on the community page that this may be because the header is missing some relevant meta-data information? Any help on how to fix this would be greatly appreciated!
This is cross-posted on the community page for ONT here:
https://community.nanoporetech.com/posts/albacore-basecalling-runni

I believe you have multiread fast5 files, as generated by recent versions of MinKNOW. A conversion script is available to convert these to the older, one-read one-file format. AFAIK Albacore is being deprecated and no version will be made supporting the multiread fast5 files, as Albacore is being replaced by Guppy.
The command line to convert from multi read fast5s to single read fast5s is as follows:
$multi_to_single_fast5 -i /path/to/fast5/ -s /path/to/save/output/

Hi Wouter, thanks for the answer! This worked as expected. To clear any ambiguity, I edited to add the command-line for future viewers. You can choose to accept or not, I just thought it added clarity! I appreciate the link.


text:
How to check if indels in VCF files are left or right aligned?
I downloaded a VCF file from dbSNP, and I'm curious if the indels in the file are left-aligned for GRCH37 genome. The documentation doesn't say anything.
How can we tell if a VCF file has left or right aligned indels?

One possibility is to run bcftools norm on the file. At the end it will print out a statistic about how many variants were realigned. I did this for dbSNPv151 (GRCh38) only for chromosome 22 like this:
$ bcftools view -i 'TYPE="indel"' dbSNP_GRCh38.p7_b151.vcf.gz 22|bcftools norm -Ou -f GRCh38.fa > /dev/null
Lines   total/split/realigned/skipped:  901503/0/490/0

Hi, although it's too late, would you please let me know which type of human reference genome, toplevel or primary_assembly or something else, is appropriate?

By definition in the VCF format specification the POS field is left aligned

POS int32 t 0-based leftmost coordinate

To confirm that you VCF is following the standard you would have to look at the reference and the match up the entries in the VCF.  However if it is from a reliable source like dbSNP you can probably assume they follow the spec


text:
data storage requirements for RNA-seq and WGS
Supposing upcoming RNA-seq for 3 conditions and 3 replications for each, how much is data storage requirements?

What sort of sequencing depth and read length?

Thank you let's say 75pb, paired end and 30X

I assume you mean 30 million reads (pairs really), since 30X is meaningless for RNA-seq.

Sorry I am very bad in these staffs :( thank you for correcting me

Please add all the extra information to your question since comments are easy to miss and can be deleted without warning. Also please clarify what species you will be working on since that can radically affect how much space you will need.

Thank you, I am not sure about how many samples and what experiments we are going to run yet. As soon as I obtained more information I will elaborate my post

The breakdown would be:

1.5-2 GB per sample for fastq files.
~3GB per BAM file, depending on the aligner

So ~5GB per BAM file, which is about 50GB total. That can vary by maybe 20%, depending on how compressible the results are (I rearrange our fastq files to make them more compressible, so I'm probably underestimating what others would see).
Obviously using Salmon/Kallisto would halve the storage requirements. Using CRAM would also decrease things, though obviously fewer tools support it.

So, for a whole genome sequencing that would roughly 15 TB?

You would need  600-750 million pairs for a human genome (mostly depending on the fragment size), which is 20-25 times as much data.

Maybe it is for a different question but how do you rearrange the fastq to make them more compressible? Is there some paper or blog to read on this?

@llrs It's a by-product of using clumpify (from BBTools) to remove optical duplicates. The process involves sorting fastq files such that similar reads are near each other, which means the resulting files compress more than they would otherwise. I run that as part of our normal Illumina demultiplexing pipeline.

Ok, then it might be out of my scope (I usually get the demultiplexed files). Thanks for the tip!


text:
scRNA-seq differential transcript usage
Many of the modern gene-quantification tools (Salmon/Kallisto) output transcript-level (as opposed to gene-level) data. All of the scRNA-seq analysis I have seen just uses the gene-level values. I suppose I can just use the transcript-level matrix, but then I lose any links between the different isoforms of the same gene. Maybe it does not really matter?
There are a few options for bulk RNA-seq and there was a similar question here, but are there any single-cell options?

Single cell data tends to be very spotty, I don't know that you'll have enough reads to do a transcript level analysis from a typical sc-RNA experiment.

@swbarnes2 True, but it's spotty even on gene level.

Could you please clarify the question? You want to know if using gene-level values in scRNA is fine with later methods or if it is fine to use transcript level values with later analysis? Or do you want to perform differential transcript usage and you don't know if it is fine (with what? other methods the conclusions the data?...) ?

@llrs I was trying to leave the question a little open-ended. The real question is at the end: what are DTU options for scRNA-seq?

I was hoping to find some tool that already exists.

The short answer is, none exist yet. I see a number of hurdles to overcome first.1
For example, no one has yet characterized saturation for transcript discovery, i.e., the reads per cell needed to maximize transcript information recovery.  This is a typical stepping-stone as a sequencing approach matures, and until it's done, it makes it unlikely that most scRNA-seq datasets are going to be appropriate vis-à-vis statistical power, even if there was a tool.  That is, if people are only aiming for cell-type discrimination via gene-level markers, datasets will likely be underpowered to measure transcript level differences.
Another issue is the lack of any transcript-level scRNA-seq simulators. All the simulators out there that I know of produce gene-level matrices.  Having a transcript-level simulator is crucial for benchmarking and characterizing the performance of differential expression calling methodologies.  Without anything to test against, that leaves a huge blindspot in terms of estimating the uncertainty (false positive rates, sensitivity) that any tool might have.  However, since some scRNA-seq models (e.g., Velten, et al. 2015, or Lopez, et al., 2018) are generative, perhaps this is shouldn't be too difficult to adapt to generating FASTQs.
Simulators are fine and all, but ground truth is probably the key hurdle of the enterprise right now. While there are definitely groups using smFISH to distinguish isoforms quantitatively on a single-cell level, there doesn't appear to be enough smFISH data out there to measure a transcriptome-wide scRNA-seq method against. Hopefully this will change in the coming year or two, as results from more highly-multiplexed single-molecule technologies become available.
1 This is just my opinion formed as I try to keep up with the literature, so I could be missing things.  Comments/corrections certainly welcome!

There is a paper that takes advantage of transcript compatibility counts (TCCs). A TCC describes all of the potential transcripts that a cDNA alignment can be derived from.  The senior author, Lior Pachter also has a nice blog post about differential transcript usage.
A discriminative learning approach to differential expression analysis for single-cell RNA-seq

Single-cell RNA-seq makes it possible to characterize the
  transcriptomes of cell types across different conditions and to
  identify their transcriptional signatures via differential analysis.
  Our method detects changes in transcript dynamics and in overall gene
  abundance in large numbers of cells to determine differential
  expression. When applied to transcript compatibility counts obtained
  via pseudoalignment, our approach provides a quantification-free
  analysis of 3′ single-cell RNA-seq that can identify previously
  undetectable marker genes.

That's a nice paper. Unfortunately, it's more of a proof-of-concept. I was hoping to find some tool that already exists.


text:
How to quantile normalization on RNA seq counts
I have a read count data (RNAseq) and want to perform quantile normalization. Could you please help me how to do it. I tried some scripts in R but it didn't work. I want the result output in matrix form.
  gene_id   SRR896664   SRR896663   SRR896665
  ENSG00000000003   46106   36353   40614
  ENSG00000000005   198 399 1200
  ENSG00000000419   40364   37769   40849
  ENSG00000000457   18924   16211   16057
  ENSG00000000460   31040   28888   29901
  ENSG00000000938   200 0   394
  ENSG00000000971   14935   14353   12522

The script which i tried is
  data <- read.csv("testquantile.csv",header=T)
  head(data)
  rownames(data) <- data[,1]
  data_mat <- data.matrix(data[,-1])
  head(data_mat)
  data_norm <- normalize.quantiles(m, copy = TRUE)

@swbarnes2 You can simply suggest to include the script, rather than giving a moral lecture. It is Bioinformatics forum, not a platform to taunt anyone.

@priya Thanks for including the script. The code seems reasonable. Did it give you an error? For what command? What was the error?

@burger Thank you for getting back. My output file was showing a result for only 2 columns.

On google there are many tutorials about quantile normalzation, for example here. In that tutorial they made a function to calculate quantile normalization. Here an example with that function on your small data set.
data
          gene_id SRR896664 SRR896663 SRR896665
1 ENSG00000000003     46106     36353     40614
2 ENSG00000000005       198       399      1200
3 ENSG00000000419     40364     37769     40849
4 ENSG00000000457     18924     16211     16057
5 ENSG00000000460     31040     28888     29901
6 ENSG00000000938       200         0       394
7 ENSG00000000971     14935     14353     12522

rownames(data) <- data$gene_id

quantile_normalisation <- function(df){
  df_rank <- apply(df,2,rank,ties.method="min")
  df_sorted <- data.frame(apply(df, 2, sort))
  df_mean <- apply(df_sorted, 1, mean)

  index_to_mean <- function(my_index, my_mean){
    return(my_mean[my_index])
  }

  df_final <- apply(df_rank, 2, index_to_mean, my_mean=df_mean)
  rownames(df_final) <- rownames(df)
  return(df_final)
}

new_data <- quantile_normalisation(data[,2:4])
new_data
                 SRR896664  SRR896663  SRR896665
ENSG00000000003 41574.6667 39110.3333 39110.3333
ENSG00000000005   197.3333   599.6667   599.6667
ENSG00000000419 39110.3333 41574.6667 41574.6667
ENSG00000000457 17064.0000 17064.0000 17064.0000
ENSG00000000460 29943.0000 29943.0000 29943.0000
ENSG00000000938   599.6667   197.3333   197.3333
ENSG00000000971 13936.6667 13936.6667 13936.6667

boxplot(data[,2:4])

boxplot(new_data)

Ma be CQN from Bioconductor will be useful, though it doesn't perform just quantile normalisation.

Welcome to the site.  Could you expand your answer to include why CQN is a good choice?

Hi, I changed my answer. I realised its not what exactly op wants.

https://www.rdocumentation.org/packages/limma/versions/3.28.14/topics/normalizeQuantiles

is a possibility.
Another popular method is https://www.rdocumentation.org/packages/preprocessCore/versions/1.34.0/topics/normalize.quantiles.
library(preprocessCore)
normalize.quantiles

the limma function you mentioned is actually used for microarray datasets not RNAseq ones.


text:
Can a gff file be converted to a fasta file?
I downloaded an annotated genome file in gff format here. I would like to use it for proteomics. Though I need it in fasta format. Is there any tool that converts gff to annotated fasta? I see the file contains the genome sequence in fasta format, but not annotated. Are those file formats convertable?

A gff is a file of annotation.  It generally doesn't include sequence information, so you can not inter-convert.  But the file on your page has a fasta entry at the bottom..  I'd just grab the bottom half of it, and if the sequence length looks right, that's likely what you want.

Only to some degree. I would like to conserve the annotations.

@Sören could you clarify what you mean by "conserve the annotations"? The annotations are in gff format, that's standard. Do you want to add the annotations to the fasta file for some reason?

You can use GFF utilities to get a fasta sequence from a GFF file using a reference genome file.

Is there any tool that converts gff to annotated fasta?

This is not an uncommon question / request, but perhaps a crash course refresher on file formats would provide some valuable context.

FASTA files are designed to hold nucleotide or amino acid sequences. The "defline" (header) of each record may include some unique identifiers or descriptive metadata. There are no strict requirements, although there are many common conventions. The information that should be the defline of a FASTA record depends primarily on the tool(s) used to produce and analyze the data.
GFF files are designed to hold annotations corresponding to sequences, most often genomic sequences such as chromosomes or (if it's a draft genome assembly) scaffolds or contigs. Each record in the GFF file specifies a sequence (such as chr6 or scaffold1985) and an interval (such as 19920557-19997344 or 2488-2489). The GFF file labels that interval as a gene or transposable element or transcription factor binding site or whatever genomic feature you're interested in.

So in this context, a "tool that converts gff to annotated fasta" doesn't make sense since it's an underdefined task. There are many ways to extract sequences from a Fasta file that has been annotated with a GFF file. Many, many scientists have written throwaway programs to do this type of task. It can be done in a hundred subtly different ways, and most of these programs are designed for a single narrow use case. Do you want to:

pull out the sequence corresponding to every GFF record? Usually not helpful, since many records are children of other records.
pull out the sequence corresponding to every gene? More common, but since you want to do proteomics you're probably not interested in promoters or introns (if either of these are annotated in your GFF file).
pull out the coding sequence of each gene and translate it into its corresponding peptide sequence?

It looks like that last point is most likely, but unless I missed something it looks like Prokka is already doing that for you. So the question remains: precisely what information exists in your GFF file that you would like in your FASTA file? You've provided an example of what you have (the FASTA files at least), but can you provide an example of what you want? Only then can we help you to identify the right tool for the task. There's a good chance it will involve writing some code in a programming language like Python or Perl.

I extracted the sequence contained in the gff file into a separate file and added a header according to the fasta file format specifications. I named that file ERS654933.fasta and then used prokka to generate my own annotation.
File header ERS654933.fasta:
>ATCC_NCTN13373
GTGTATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGG
GAAAAAGTGCTTGAAATTGCTCAAGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTA
AAAGATACTGAGCTTTACACGATCAAAGATGGTGAAGCTATCGTATTATCGAGTATTCCT
TTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGAT
GTTGTAGGCTATGAAGTAAAACCTCACTTTATTACTACTGAAGAATTAGCAAATTATAGT

Then:
prokka ERS654933.fasta

Which generated a folder with several files. Here is the header of the generated file PROKKA_01232019.faa.
NEOBMGJD_00001 Chromosomal replication initiator protein DnaA
MSEKEIWEKVLEIAQEKLSAVSYSTFLKDTELYTIKDGEAIVLSSIPFNANWLNQQYAEI
IQAILFDVVGYEVKPHFITTEELANYSNNETATPKEATKPSTETTEDNHVLGREQFNAHN
TFDTFVIGPGNRFPHAASLAVAEAPAKAYNPLFIYGGVGLGKTHLMHAIGHHVLDNNPDA
KVIYTSSEKFTNEFIKSIRDNEGEAFRERYRNIDVLLIDDIQFIQNKVQTQEEFFYTFNE
LHQNNKQIVISSDRPPKEIAQLEDRLRSRFEWGLIVDITPPDYETRMAILQKKIEEEKLD
IPPEALNYIANQIQSNIRELEGALTRLLAYSQLLGKPITTELTAEALKDIIQAPKSKKIT
IQDIQKIVGQYYNVRIEDFSAKKRTKSIAYPRQIAMYLSRELTDFSLPKIGEEFGGRDHT
TVIHAHEKISKDLKEDPIFKQEVENLEKEIRNV
>NEOBMGJD_00002 Beta sliding clamp
MMEFTIKRDYFITQLNDTLKAISPRTTLPILTGIKIDAKEHEVILTGSDSEISIEITIPK
TVDGEDIVNISETGSVVLPGRFFVDIIKKLPGKDVKLSTNEQFQTLITSGHSEFNLSGLD
PDQYPLLPQVSRDDAIQLSVKVLKNVIAQTNFAVSTSETRPVLTGVNWLIQENELICTAT
DSHRLAVRKLQLEDVSENKNVIIPGKALAELNKIMSDNEEDIDIFFASNQVLFKVGNVNF
ISRLLEGHYPDTTRLFPENYEIKLSIDNGEFYHAIDRASLLAREGGNNVIKLSTGDDVVE
LSSTSPEIGTVKEEVDANDVEGGSLKISFNSKYMMDALKAIDNDEVEVEFFGTMKPFILK
PKGDDSVTQLILPIRTY
>NEOBMGJD_00003 hypothetical protein
MVQEVVVEGDINLGQFLKTEGIIESGGQAKWFLQDVEVLINGVRETRRGKKLEHQDRIDI
PELPEDAGSFLIIHQGEQ

So, that way I could create my own annotation, but I am loosing the annotation information in the gff file.

You've got gene name and description, and you want more?  You might have to make your own script for that.

I thought, I am probably not the first one doing that. So, I was hoping someone knows how to do that. Rather then inventing the wheel again. I have done that before.