# Question
Title: Duplicate genes with RSEM counts: Which one to choose?
I have Ensembl ids in the first column and samples with RSEM counts data in other columns. I converted Ensembl ids to gene symbols. Now I see there are three genes repeated twice.

```
Gene   S1       S2       S3
COG8  804.07  1475.16  323.80
COG8  249.23   867.94  113.64
SDHD 1464.84  3508.02    75
SDHD  4.16     3.97       0
SOGA3 81.65   30.76    32.69
SOGA3 0.00    10.77    25.65

```

These are the Ensembl ids:

* SOGA3: ENSG00000214338, ENSG00000255330
* SDHD: ENSG00000204370, ENSG00000255292
* COG8: ENSG00000213380, ENSG00000272617.

This is not for DE analysis, I want to do subtyping. Should I calculate mean or median on that? How should I go forward with this duplicate genes?

# Answer
> 5 votes
Let's look into this a bit more deeply. For instance:

The Ensembl pages (linked above) for both ENSG00000214338 and ENSG00000255330 show that these genes have 3 paralogs and each lists the other as one of the paralogs. Your other two examples are similar, they are also paralogs of each other. You can also see that the position of each gene in the genome is slightly different:

ENSG00000255330 : Chromosome 6: 127,438,406-127,519,001 reverse strand. ENSG00000214338 : Chromosome 6: 127,472,794-127,519,191 reverse strand.

So these two genes refer to two distinct regions of the genome. And, as expected, they have different transcripts. However, two of those transcripts (ENST00000481848.6 from ENSG00000255330 and ENST00000525778.5 from ENSG00000214338) code for the same protein sequence. This means that we have *one* protein but *two* different genes for it. To illustrate, here are the two transcripts aligned:

```
ENST00000525778.5       CCGTGTTCGAGCGGTTTTGGGATCCTCTGCCAAGTACGGGATAAGAGAAT
ENST00000481848.6       CTC-----------------------------------------------
                        *

ENST00000525778.5       AGAGCAAGGAGTTTGAGCGCTGCTGCCAATCTTCCATCTCGGGCGTGGCT
ENST00000481848.6       --------------------------------------------------

ENST00000525778.5       CTGGCCTTTTTGTCTCTTTATCCCGCCACTCCCCACCCCCGTCCCCCGTC
ENST00000481848.6       --------------------------------------------------

ENST00000525778.5       CCCCTCCCCGCCTTCTTCTCTTTCCTCTTGAAAGTAACCTCTCTTTGGTT
ENST00000481848.6       -------------------------------------------TTTGGTT
                                                                   *******

ENST00000525778.5       ACTTCTACCCTATCGTCCCTTGCAATCCCGCCTTCTCTCCTGCCCAGTTC
ENST00000481848.6       ACTTCTACCCTATCGTCCCTTGCAATCCCGCCTTCTCTCCTGCCCAGTTC
                        **************************************************

ENST00000525778.5       GGTGAGGTCAGCTCATGAACATCTCTGCACTTCCTGGGGAGCGTCTGGCC
ENST00000481848.6       GGTGAGGTCAGCTCATGAACATCTCTGCACTTCCTGGGGAGCGTCTGGCC
                        **************************************************

ENST00000525778.5       CAGCTTGGTACCTCGCGTGTTAGCCTGCTGAGAGCTGCAGCTCCAGCCCA
ENST00000481848.6       CAGCTTGGTACCTCGCGTGTTAGCCTGCTGAGAGCTGCAGCTCCAGCCCA
                        **************************************************

ENST00000525778.5       GCCGCGGACCGCCGGGAAGTCACCCGCTCTGGGGGCAGCTGAATCTTGAC
ENST00000481848.6       GCCGCGGACCGCCGGGAAGTCACCCGCTCTGGGGGCAGCTGAATCTTGAC
                        **************************************************

ENST00000525778.5       GTTACCGCCATTGTTTGGAGAAATCGAGTGGGTGTGGAACTGCGGCGAGC
ENST00000481848.6       GTTACCGCCATTGTTTGGAGAAATCGAGTGGGTGTGGAACTGCGGCGAGC
                        **************************************************

ENST00000525778.5       GAGCCCAAAGTCATGCTCAGCCGGAAAGTCGCATACAATAATATGAAGTT
ENST00000481848.6       GAGCCCAAAGTCATGCTCAGCCGGAAAGTCGCATACAATAATATGAAGTT
                        **************************************************

ENST00000525778.5       GAACTTTCCTCGGCGCAGAAAGGCACAATAATAAACCCTGAGATCCCCGG
ENST00000481848.6       GAACTTTCCTCGGCGCAGAAAGGCACAATAATAAACCCTGAGATCCCCGG
                        **************************************************

ENST00000525778.5       AATATCAGAAGGAAAGTGCTGTTGATTTGTTCAGCTACTTGGAAACCTCG
ENST00000481848.6       AATATCAGAAGGAAAGTGCTGTTGATTTGTTCAGCTACTTGGAAACCTCG
                        **************************************************

ENST00000525778.5       CCTCTGCTGAGAAGTGAGCAAGTTCTGGATTTCTAAAGGGTTGGGGATGA
ENST00000481848.6       CCTCTGCTGAGAA-------------------------------------
                        *************

ENST00000525778.5       ATGATGGGCTTAGGAAGGAAGGTTGGCGGCGGTGAAGAGGGAAGAGGATG
ENST00000481848.6       --------CTTAGGAAGGAAGGTTGGCGGCGGTGAAGAGGGAAGAGGATG
                                ******************************************

ENST00000525778.5       ATCTAGTAATGGGCTCTTCCAAATCGCATCTGGTAGGTTCATAGCCATGA
ENST00000481848.6       ATCTAGTAATGGGCTCTTCCAAATCGCATCTGGTAGGTTCATAGCCATGA
                        **************************************************

ENST00000525778.5       GTCAGCCACCAATCGGGGGCGCTGCCCCTGCCACAGCAGCAGCTTCTCCC
ENST00000481848.6       GTCAGCCACCAATCGGGGGCGCTGCCCCTGCCACAGCAGCAGCTTCTCCC
                        **************************************************

ENST00000525778.5       GCCGCTGCAGCCACTGAGGCTCGTCTGCACCCGGAGGGCAGCAGCAGAAA
ENST00000481848.6       GCCGCTGCAGCCACTGAGGCTCGTCTGCACCCGGAGGGCAGCAGCAGAAA
                        **************************************************

ENST00000525778.5       GCAGCAGCGAGCTCAGTCACCCGCTAGACCGAGGGACAGTTCGCTGCGAC
ENST00000481848.6       GCAGCAGCGAGCTCAGTCACCCGCTAGACCGAGGGACAGTTCGCTGCGAC
                        **************************************************

ENST00000525778.5       AGACAATCGCGGCCACCCGGTCCCCAGTGGGGGCCGGCACTAAACTCAAT
ENST00000481848.6       AGACAATCGCGGCCACCCGGTCCCCAGTGGGGGCCGGCACTAAACTCAAT
                        **************************************************

ENST00000525778.5       TCTGTTCGGCAGCAGCAGCTGCAGCAGCAGCAACAGCAGGGTAACAAGAC
ENST00000481848.6       TCTGTTCGGCAGCAGCAGCTGCAGCAGCAGCAACAGCAGGGTAACAAGAC
                        **************************************************

ENST00000525778.5       CGGGAGCCGCACGGGGCCGCCCGCCAGCATCCGAGGAGGCGGAGGTGGGG
ENST00000481848.6       CGGGAGCCGCACGGGGCCGCCCGCCAGCATCCGAGGAGGCGGAGGTGGGG
                        **************************************************

ENST00000525778.5       CCGAGAAGGCGACGCCCCTGGCGCCCAAAGGGGCCGCTCCGGGAGCTGTC
ENST00000481848.6       CCGAGAAGGCGACGCCCCTGGCGCCCAAAGGGGCCGCTCCGGGAGCTGTC
                        **************************************************

ENST00000525778.5       CAGCCTGTGGCTGGTGCTGAAGCGGCCCCCGCCGCGACCCTGGCCGCTTT
ENST00000481848.6       CAGCCTGTGGCTGGTGCTGAAGCGGCCCCCGCCGCGACCCTGGCCGCTTT
                        **************************************************

ENST00000525778.5       GGGCGGCAGGAGGCCCGGACCGCCCGAGGAGCCGCCTCGGGAGCTAGAGT
ENST00000481848.6       GGGCGGCAGGAGGCCCGGACCGCCCGAGGAGCCGCCTCGGGAGCTAGAGT
                        **************************************************

ENST00000525778.5       CCGTGCCCTCGAAACTCGGGGAACCCCCTCCGCTCGGAGAAGGGGGAGGA
ENST00000481848.6       CCGTGCCCTCGAAACTCGGGGAACCCCCTCCGCTCGGAGAAGGGGGAGGA
                        **************************************************

ENST00000525778.5       GGGGGCGGGGAAGGAGGAGGAGCCGGCGGCGGCTCCGGGGAAAGGGAGGG
ENST00000481848.6       GGGGGCGGGGAAGGAGGAGGAGCCGGCGGCGGCTCCGGGGAAAGGGAGGG
                        **************************************************

ENST00000525778.5       GGGCGCTCCGCAGCCGCCGCCGCCCAGGGGCTGGCGAGGGAAAGGCGTAC
ENST00000481848.6       GGGCGCTCCGCAGCCGCCGCCGCCCAGGGGCTGGCGAGGGAAAGGCGTAC
                        **************************************************

ENST00000525778.5       GCGCTCAGCAGAGGGGCGGCAGCGGCGGGGAGGGGGCCTCCCCTTCTCCA
ENST00000481848.6       GCGCTCAGCAGAGGGGCGGCAGCGGCGGGGAGGGGGCCTCCCCTTCTCCA
                        **************************************************

ENST00000525778.5       TCCTCCTCTTCTGCGGGCAAAACCCCAGGAACCGGCAGCAGAAACTCCGG
ENST00000481848.6       TCCTCCTCTTCTGCGGGCAAAACCCCAGGAACCGGCAGCAGAAACTCCGG
                        **************************************************

ENST00000525778.5       AAGCGGCGTTGCGGGGGGCGGCAGCGGTGGTGGAGGGAGCTACTGGAAAG
ENST00000481848.6       AAGCGGCGTTGCGGGGGGCGGCAGCGGTGGTGGAGGGAGCTACTGGAAAG
                        **************************************************

ENST00000525778.5       AAGGATGTCTGCAGTCTGAGCTCATCCAGTTCCATCTCAAGAAGGAGCGG
ENST00000481848.6       AAGGATGTCTGCAGTCTGAGCTCATCCAGTTCCATCTCAAGAAGGAGCGG
                        **************************************************

ENST00000525778.5       GCGGCAGCGGCGGCGGCCGCGGCTCAGATGCACGCTAAGAACGGCGGCGG
ENST00000481848.6       GCGGCAGCGGCGGCGGCCGCGGCTCAGATGCACGCTAAGAACGGCGGCGG
                        **************************************************

ENST00000525778.5       CAGCAGTAGCCGCAGCTCCCCGGTGTCTGGCCCCCCTGCCGTTTGCGAGA
ENST00000481848.6       CAGCAGTAGCCGCAGCTCCCCGGTGTCTGGCCCCCCTGCCGTTTGCGAGA
                        **************************************************

ENST00000525778.5       CCCTGGCCGTCGCCTCCGCCTCCCCAATGGCGGCGGCGGCGGAGGGCCCC
ENST00000481848.6       CCCTGGCCGTCGCCTCCGCCTCCCCAATGGCGGCGGCGGCGGAGGGCCCC
                        **************************************************

ENST00000525778.5       CAGCAGAGCGCAGAGGGCAGCGCGAGCGGCGGGGGCATGCAGGCGGCAGC
ENST00000481848.6       CAGCAGAGCGCAGAGGGCAGCGCGAGCGGCGGGGGCATGCAGGCGGCAGC
                        **************************************************

ENST00000525778.5       GCCCCCTTCGTCGCAGCCGCACCCGCAGCAGCTCCAAGAGCAGGAAGAAA
ENST00000481848.6       GCCCCCTTCGTCGCAGCCGCACCCGCAGCAGCTCCAAGAGCAGGAAGAAA
                        **************************************************

ENST00000525778.5       TGCAAGAGGAGATGGAGAAGCTGCGAGAGGAAAACGAGACTCTCAAGAAC
ENST00000481848.6       TGCAAGAGGAGATGGAGAAGCTGCGAGAGGAAAACGAGACTCTCAAGAAC
                        **************************************************

ENST00000525778.5       GAGATCGATGAGCTGAGAACCGAGATGGACGAGATGAGGGACACTTTCTT
ENST00000481848.6       GAGATCGATGAGCTGAGAACCGAGATGGACGAGATGAGGGACACTTTCTT
                        **************************************************

ENST00000525778.5       CGAGGAGGATGCCTGTCAACTGCAGGAAATGCGCCACGAGTTGGAGAGAG
ENST00000481848.6       CGAGGAGGATGCCTGTCAACTGCAGGAAATGCGCCACGAGTTGGAGAGAG
                        **************************************************

ENST00000525778.5       CCAACAAAAACTGCCGGATCCTGCAGTACCGCCTCCGCAAAGCCGAGCGC
ENST00000481848.6       CCAACAAAAACTGCCGGATCCTGCAGTACCGCCTCCGCAAAGCCGAGCGC
                        **************************************************

ENST00000525778.5       AAAAGGCTCCGCTACGCCCAGACCGGGGAAATCGACGGGGAGCTGTTGCG
ENST00000481848.6       AAAAGGCTCCGCTACGCCCAGACCGGGGAAATCGACGGGGAGCTGTTGCG
                        **************************************************

ENST00000525778.5       CAGCCTGGAGCAGGACCTCAAGGTTGCAAAGGATGTATCTGTGAGACTTC
ENST00000481848.6       CAGCCTGGAGCAGGACCTCAAGGTTGCAAAGGATGTATCTGTGAGACTTC
                        **************************************************

ENST00000525778.5       ACCATGAATTAGAAAATGTGGAAGAAAAGAGAACAACAACAGAAGATGAA
ENST00000481848.6       ACCATGAATTAGAAAATGTGGAAGAAAAGAGAACAACAACAGAAGATGAA
                        **************************************************

ENST00000525778.5       AATGAGAAACTGAGGCAACAGCTCATAGAAGTTGAAATTGCAAAGCAAGC
ENST00000481848.6       AATGAGAAACTGAGGCAACAGCTCATAGAAGTTGAAATTGCAAAGCAAGC
                        **************************************************

ENST00000525778.5       TTTACAGAATGAACTGGAAAAAATGAAAGAGTTATCCTTAAAAAGAAGAG
ENST00000481848.6       TTTACAGAATGAACTGGAAAAAATGAAAGAGTTATCCTTAAAAAGAAGAG
                        **************************************************

ENST00000525778.5       GAAGCAAAGATTTGCCAAAATCTGAAAAAAAGGCTCAACAGACTCCCACA
ENST00000481848.6       GAAGCAAAGATTTGCCAAAATCTGAAAAAAAGGCTCAACAGACTCCCACA
                        **************************************************

ENST00000525778.5       GAGGAGGACAATGAAGATCTGAAGTGCCAGCTGCAGTTTGTTAAGGAAGA
ENST00000481848.6       GAGGAGGACAATGAAGATCTGAAGTGCCAGCTGCAGTTTGTTAAGGAAGA
                        **************************************************

ENST00000525778.5       AGCCGCTTTGATGAGAAAGAAAATGGCCAAGATTGATAAAGAAAAGGACA
ENST00000481848.6       AGCCGCTTTGATGAGAAAGAAAATGGCCAAGATTGATAAAGAAAAGGACA
                        **************************************************

ENST00000525778.5       GATTTGAACACGAGCTCCAGAAGTACAGATCCTTTTATGGGGATCTGGAC
ENST00000481848.6       GATTTGAACACGAGCTCCAGAAGTACAGATCCTTTTATGGGGATCTGGAC
                        **************************************************

ENST00000525778.5       AGTCCTTTGCCCAAAGGAGAAGCCGGAGGCCCTCCCAGCACTAGGGAGGC
ENST00000481848.6       AGTCCTTTGCCCAAAGGAGAAGCCGGAGGCCCTCCCAGCACTAGGGAGGC
                        **************************************************

ENST00000525778.5       CGAGCTCAAGCTACGGCTAAGGCTGGTGGAGGAAGAAGCCAACATCCTGG
ENST00000481848.6       CGAGCTCAAGCTACGGCTAAGGCTGGTGGAGGAAGAAGCCAACATCCTGG
                        **************************************************

ENST00000525778.5       GCAGGAAAATCGTCGAACTGGAGGTGGAGAACAGAGGCCTGAAGGCGGAA
ENST00000481848.6       GCAGGAAAATCGTCGAACTGGAGGTGGAGAACAGAGGCCTGAAGGCGGAA
                        **************************************************

ENST00000525778.5       CTGGACGACCTTAGGGGCGATGACTTCAACGGCTCGGCCAACCCGCTCAT
ENST00000481848.6       CTGGACGACCTTAGGGGCGATGACTTCAACGGCTCGGCCAACCCGCTCAT
                        **************************************************

ENST00000525778.5       GAGGGAGCAGAGCGAATCCCTGTCGGAGCTGCGGCAGCACCTGCAGCTGG
ENST00000481848.6       GAGGGAGCAGAGCGAATCCCTGTCGGAGCTGCGGCAGCACCTGCAGCTGG
                        **************************************************

ENST00000525778.5       TGGAAGACGAGACGGAGCTGCTGCGGAGGAACGTGGCCGACCTGGAGGAG
ENST00000481848.6       TGGAAGACGAGACGGAGCTGCTGCGGAGGAACGTGGCCGACCTGGAGGAG
                        **************************************************

ENST00000525778.5       CAGAACAAGCGCATCACGGCGGAGCTCAACAAGTACAAGTACAAGTCCGG
ENST00000481848.6       CAGAACAAGCGCATCACGGCGGAGCTCAACAAGTACAAGTACAAGTCCGG
                        **************************************************

ENST00000525778.5       CGGCCACGACAGCGCGCGGCACCACGACAACGCCAAGACCGAGGCCCTGC
ENST00000481848.6       CGGCCACGACAGCGCGCGGCACCACGACAACGCCAAGACCGAGGCCCTGC
                        **************************************************

ENST00000525778.5       AGGAGGAGCTGAAGGCGGCGCGCCTGCAGATCAACGAGCTCAGCGGCAAG
ENST00000481848.6       AGGAGGAGCTGAAGGCGGCGCGCCTGCAGATCAACGAGCTCAGCGGCAAG
                        **************************************************

ENST00000525778.5       GTCATGCAGCTGCAGTACGAGAACCGCGTGCTTATGTCCAACATGCAGCG
ENST00000481848.6       GTCATGCAGCTGCAGTACGAGAACCGCGTGCTTATGTCCAACATGCAGCG
                        **************************************************

ENST00000525778.5       CTACGACCTGGCCTCGCACCTGGGCATCCGCGGCAGCCCCCGCGACAGCG
ENST00000481848.6       CTACGACCTGGCCTCGCACCTGGGCATCCGCGGCAGCCCCCGCGACAGCG
                        **************************************************

ENST00000525778.5       ACGCCGAGAGCGACGCGGGCAAGAAGGAGAGCGACGACGACTCGCGGCCT
ENST00000481848.6       ACGCCGAGAGCGACGCGGGCAAGAAGGAGAGCGACGACGACTCGCGGCCT
                        **************************************************

ENST00000525778.5       CCGCACCGCAAGCGCGAAGGGCCCATCGGCGGCGAGAGCGACTCGGAGGA
ENST00000481848.6       CCGCACCGCAAGCGCGAAGGGCCCATCGGCGGCGAGAGCGACTCGGAGGA
                        **************************************************

ENST00000525778.5       GGTGCGCAACATCCGCTGCCTCACGCCCACTCGCTCCTTCTACCCGGCGC
ENST00000481848.6       GGTGCGCAACATCCGCTGCCTCACGCCCACTCGCTCCTTCTACCCGGCGC
                        **************************************************

ENST00000525778.5       CCGGGCCCTGGCCCAAGAGCTTCTCCGATCGGCAGCAGATGAAGGACATC
ENST00000481848.6       CCGGGCCCTGGCCCAAGAGCTTCTCCGATCGGCAGCAGATGAAGGACATC
                        **************************************************

ENST00000525778.5       CGCTCGGAGGCCGAGCGCCTGGGCAAGACCATCGACCGGCTCATCGCCGA
ENST00000481848.6       CGCTCGGAGGCCGAGCGCCTGGGCAAGACCATCGACCGGCTCATCGCCGA
                        **************************************************

ENST00000525778.5       CACGAGCACCATCATCACCGAGGCGCGCATCTACGTGGCCAACGGGGACC
ENST00000481848.6       CACGAGCACCATCATCACCGAGGCGCGCATCTACGTGGCCAACGGGGACC
                        **************************************************

ENST00000525778.5       TGTTCGGACTCATGGACGAGGAGGACGACGGCAGCCGCATCCGGGAGCAC
ENST00000481848.6       TGTTCGGACTCATGGACGAGGAGGACGACGGCAGCCGCATCCGGGAGCAC
                        **************************************************

ENST00000525778.5       GAGCTGCTCTACCGCATCAACGCTCAGATGAAGGCCTTCCGCAAGGAGCT
ENST00000481848.6       GAGCTGCTCTACCGCATCAACGCTCAGATGAAGGCCTTCCGCAAGGAGCT
                        **************************************************

ENST00000525778.5       GCAGACCTTCATCGACCGCCTCGAGGTGCCCAAGTCTGCGGACGACCGCG
ENST00000481848.6       GCAGACCTTCATCGACCGCCTCGAGGTGCCCAAGTCTGCGGACGACCGCG
                        **************************************************

ENST00000525778.5       GCGCCGAGGAGCCCATTTCCGTGAGTCAGATGTTCCAGCCTATCATTTTA
ENST00000481848.6       GCGCCGAGGAGCCCATTTCCGTGAGTCAGATGTTCCAGCCTATCATTTTA
                        **************************************************

ENST00000525778.5       CTTATTCTCATTCTTGTATTATTTTCATCACTTTCTTACACAACAATATT
ENST00000481848.6       CTTATTCTCATTCTTGTATTATTTTCATCACTTTCTTACACAACAATATT
                        **************************************************

ENST00000525778.5       TAAACTTGTCTTCCTTTTTACACTGTTTTTTGTACTGTAAATCTTTCATC
ENST00000481848.6       TAAACTTGTCTTCCTTTTTACACTGTTTTTTGTACTGTAAATCTTTCATC
                        **************************************************

ENST00000525778.5       ATTTACCATTCATTGTAGTATTTTCAGTTTGTTTATTTTGTTCACCCTTC
ENST00000481848.6       ATTTACCATTCATTGTAGTATTTTCAGTTTGTTTATTTTGTTCACCCTTC
                        **************************************************

ENST00000525778.5       AAGACAAGAAGTAAAAGAAGTATAATTTCTGTAGTAACCAATGCTATAAA
ENST00000481848.6       AAGACAAGAAGTAAAAGAAGTATAATTTCTGTAGTAACCAATGCTATAAA
                        **************************************************

ENST00000525778.5       AACACTGAAGACTGCTTATTTCTTTAAAAAGATACAACTCATCTTACCAA
ENST00000481848.6       AACACTGAAGACTGCTTATTTCTTTAAAAAGATACAACTCATCTTACCAA
                        **************************************************

ENST00000525778.5       GACCAAA-------------------------------------------
ENST00000481848.6       GACCAAATTCAATAAGAAGCCCAAACACTAAAATATTTCAGGCTTTATTT
                        *******

ENST00000525778.5       --------------------------------------------------
ENST00000481848.6       TAAAGGCAAGTGAGACTGCTTCAAATAAAACAACTTCAAGCTTCCAAGAA

[. . .]

```

<sub>The ENST00000481848.6 transcript continues quite a bit after this but there is no more similarity. </sub>

The transcripts are not identical, but are very similar. Now, look at the proteins each transcript encodes:

```
ENST00000481848.6  MSQPPIGGAAPATAAASPAAAATEARLHPEGSSRKQQRAQSPARPRDSSL
ENST00000525778.5  MSQPPIGGAAPATAAASPAAAATEARLHPEGSSRKQQRAQSPARPRDSSL
                   **************************************************

ENST00000481848.6  RQTIAATRSPVGAGTKLNSVRQQQLQQQQQQGNKTGSRTGPPASIRGGGG
ENST00000525778.5  RQTIAATRSPVGAGTKLNSVRQQQLQQQQQQGNKTGSRTGPPASIRGGGG
                   **************************************************

ENST00000481848.6  GAEKATPLAPKGAAPGAVQPVAGAEAAPAATLAALGGRRPGPPEEPPREL
ENST00000525778.5  GAEKATPLAPKGAAPGAVQPVAGAEAAPAATLAALGGRRPGPPEEPPREL
                   **************************************************

ENST00000481848.6  ESVPSKLGEPPPLGEGGGGGGEGGGAGGGSGEREGGAPQPPPPRGWRGKG
ENST00000525778.5  ESVPSKLGEPPPLGEGGGGGGEGGGAGGGSGEREGGAPQPPPPRGWRGKG
                   **************************************************

ENST00000481848.6  VRAQQRGGSGGEGASPSPSSSSAGKTPGTGSRNSGSGVAGGGSGGGGSYW
ENST00000525778.5  VRAQQRGGSGGEGASPSPSSSSAGKTPGTGSRNSGSGVAGGGSGGGGSYW
                   **************************************************

ENST00000481848.6  KEGCLQSELIQFHLKKERAAAAAAAAQMHAKNGGGSSSRSSPVSGPPAVC
ENST00000525778.5  KEGCLQSELIQFHLKKERAAAAAAAAQMHAKNGGGSSSRSSPVSGPPAVC
                   **************************************************

ENST00000481848.6  ETLAVASASPMAAAAEGPQQSAEGSASGGGMQAAAPPSSQPHPQQLQEQE
ENST00000525778.5  ETLAVASASPMAAAAEGPQQSAEGSASGGGMQAAAPPSSQPHPQQLQEQE
                   **************************************************

ENST00000481848.6  EMQEEMEKLREENETLKNEIDELRTEMDEMRDTFFEEDACQLQEMRHELE
ENST00000525778.5  EMQEEMEKLREENETLKNEIDELRTEMDEMRDTFFEEDACQLQEMRHELE
                   **************************************************

ENST00000481848.6  RANKNCRILQYRLRKAERKRLRYAQTGEIDGELLRSLEQDLKVAKDVSVR
ENST00000525778.5  RANKNCRILQYRLRKAERKRLRYAQTGEIDGELLRSLEQDLKVAKDVSVR
                   **************************************************

ENST00000481848.6  LHHELENVEEKRTTTEDENEKLRQQLIEVEIAKQALQNELEKMKELSLKR
ENST00000525778.5  LHHELENVEEKRTTTEDENEKLRQQLIEVEIAKQALQNELEKMKELSLKR
                   **************************************************

ENST00000481848.6  RGSKDLPKSEKKAQQTPTEEDNEDLKCQLQFVKEEAALMRKKMAKIDKEK
ENST00000525778.5  RGSKDLPKSEKKAQQTPTEEDNEDLKCQLQFVKEEAALMRKKMAKIDKEK
                   **************************************************

ENST00000481848.6  DRFEHELQKYRSFYGDLDSPLPKGEAGGPPSTREAELKLRLRLVEEEANI
ENST00000525778.5  DRFEHELQKYRSFYGDLDSPLPKGEAGGPPSTREAELKLRLRLVEEEANI
                   **************************************************

ENST00000481848.6  LGRKIVELEVENRGLKAELDDLRGDDFNGSANPLMREQSESLSELRQHLQ
ENST00000525778.5  LGRKIVELEVENRGLKAELDDLRGDDFNGSANPLMREQSESLSELRQHLQ
                   **************************************************

ENST00000481848.6  LVEDETELLRRNVADLEEQNKRITAELNKYKYKSGGHDSARHHDNAKTEA
ENST00000525778.5  LVEDETELLRRNVADLEEQNKRITAELNKYKYKSGGHDSARHHDNAKTEA
                   **************************************************

ENST00000481848.6  LQEELKAARLQINELSGKVMQLQYENRVLMSNMQRYDLASHLGIRGSPRD
ENST00000525778.5  LQEELKAARLQINELSGKVMQLQYENRVLMSNMQRYDLASHLGIRGSPRD
                   **************************************************

ENST00000481848.6  SDAESDAGKKESDDDSRPPHRKREGPIGGESDSEEVRNIRCLTPTRSFYP
ENST00000525778.5  SDAESDAGKKESDDDSRPPHRKREGPIGGESDSEEVRNIRCLTPTRSFYP
                   **************************************************

ENST00000481848.6  APGPWPKSFSDRQQMKDIRSEAERLGKTIDRLIADTSTIITEARIYVANG
ENST00000525778.5  APGPWPKSFSDRQQMKDIRSEAERLGKTIDRLIADTSTIITEARIYVANG
                   **************************************************

ENST00000481848.6  DLFGLMDEEDDGSRIREHELLYRINAQMKAFRKELQTFIDRLEVPKSADD
ENST00000525778.5  DLFGLMDEEDDGSRIREHELLYRINAQMKAFRKELQTFIDRLEVPKSADD
                   **************************************************

ENST00000481848.6  RGAEEPISVSQMFQPIILLILILVLFSSLSYTTIFKLVFLFTLFFVL
ENST00000525778.5  RGAEEPISVSQMFQPIILLILILVLFSSLSYTTIFKLVFLFTLFFVL
                   ***********************************************

```

As you can see, they are 100% identical. So, we have two genes (distinct regions of the DNA) which produce the same protein. Because the gene product is the same, the "function" of the genes is also the same so HUGO treats them as the same gene and gives them the same name.

Whether or not you want to treat them as one gene or two will depend on what you are doing. If you are comparing the relative abundance of each ***transcript***, you must treat them separately. If you want to examine the expression level of the ***protein product*** and are using the RNASEM data as a proxy for protein information, then you should treat them as one "gene" and simply sum your RNASEM values.

# Answer
> 5 votes
There is no one-to-one mapping of gene ids from one database to the other. Ensembl (who maintain Ensembl ENSG IDs), ncbi (who maintain EntrezGene IDs and RefSeq transcript ids) and HUGO (who maintain gene symbols/names) all have different ideas about what a gene is, which part of the genome belongs to which gene etc.

In fact, it's worse than that because HUGO gene IDs refer to the abstract idea of a gene, rather than a particular part of any chromosome.

My suggestion would be to leave the IDs as Ensembl IDs (as these actually refer to a physical part of a chromosome, which is what you are quantifying), right up until you have to report the results to someone who really NEEDS gene names. If you still have duplicates, you could sum, or pick the most interesting, but I would suggest just leaving the duplicates in there, and treat it as thought there are two genes with the same name.

---
Tags: r, rna-seq, bioconductor, rsem
---


# Question
Title: Using pythoncyc to get ec-code of a reaction
I am trying to collect the ec-code (and gene) of a reaction. I have tried the following:

```
import pythoncyc
pgdb = pythoncyc.PGDB("sco")
r = "RXN1A0-6128"
gene = "SCP1.233B"
enzyme = pgdb.enzymes_of_reaction(r)
enzyme = pgdb.enzymes_of_gene(gene)
pgdb.full_enzyme_name(enzyme)

```

However, I only get an enzyme id, and not the ec-code, which is visible on biocyc.

# Answer
BioGeek is correct, however I will post a full answer.

```
import pythoncyc
# To get ec-code
db = pythoncyc.select_organism('sco')
reaction = db["RXN1A0-6128"]
ec_code = reaction["ec_number"]

# To get genes
pgdb = pythoncyc.PGDB("sco")
genes = pgdb.genes_of_reaction("RXN1A0-6128")

```
> 1 votes

---
Tags: python
---


title: Proteins that cannot form biofilm?

text:
# Question
Title: Proteins that cannot form biofilm?
I am trying to build a machine learning training set for bacterial protein sequences that form biofilm, and that cannot. I collected the positive sequences from the GO ontology website but for negative sequences I am not sure which sequences to incorporate into my training set.

Can anyone point me to resources for proteins sequences that are known to be not capable of forming biofilms?

# Answer
> 2 votes
One option is to treat all proteins that are not explicitly known to form a biofilm as proteins that are incapable of doing so, but this will likely result in many false negatives.

You could take a look at one-class classification and positive-unlabeled (PU) learning. Both are techniques specifically designed for these kind of problems where only one of the two classes is known. (In this case: only biofilm forming proteins are known.)

---
Tags: proteins, machine-learning
---


# Question
Title: Estimate the length of poly-A tails from randomly-primed RNAseq data
So a poly-A tail is a long chain of adenine nucleotides that is added to a messenger RNA (mRNA) molecule during RNA processing to increase the stability of the molecule.

For my project, I would like to estimate the length of poly-A tails from randomly-primed RNAseq data.

Any ideas how this could be done? Are there algorithms/tools which currently do this?

Perhaps I would need to model this somehow? All suggestions appreciated

# Answer
> 3 votes
Just look for polyA tracts at the end of sequences, and count them if they're larger than ~18 bases. I've done this with MinION cDNA reads by mapping the polyA adapter sequence (with an elongated polyA sequence) to the reads using LAST, and working out the length of aligned sequence.

Unfortunately, because the cDNA adapter finishes with a TVN sequence, rather than TTT, binding of the adapter is not necessarily a good representation of the true polyA length. I expect that Illumina adapters will have similar problems, if not worse (due to the read length limits).

It might be possible to use direct-RNA sequencing to better categorise the polyA lengths, but I've yet to look at the public direct-RNA datasets. Just in case you're interested in looking yourself, here's a direct RNA dataset for NA12878 from the wonderful nanopore-WGS consortium:

https://github.com/nanopore-wgs-consortium/NA12878/blob/master/RNA.md

Here's the consortium's graph from that summary document showing polyA length:

# Answer
> 4 votes
Can’t be done. If you already sequenced then I’m afraid the money is wasted (unless of course the data is good for something else).

The standard Illumina basecaller doesn’t deal well with homopolymers and will consequently call incorrect poly(A) lengths.

Amongst others, Narry Kim’s lab developed a dedicated method — TAIL-seq — to solve this problem. It’s a protocol for library prep using pulldown of biotin tagged tails (rather than random priming), plus dedicated basecalling software that uses spike-ins for calibration.

Be warned though: the protocol is reportedly technically challenging, and in order to use it you need an Illumina machine that gives you access to the raw image files *before basecalling*. Up-to-date Illumina machines won’t give you these files, full stop. So you’ll need to find a machine that runs an outdated firmware.

---
Tags: rna-seq, software-recommendation, rna, algorithms
---


# Question
Title: Comparison of gene set enrichment statistics
I am performing a gene set enrichment analysis to determine if particular gene sets are coherently up- or down-regulated. I have seen several statistics for computing a p-value of GSEA-style enrichment. I'm particularly interested in the differences between a mean-rank gene set test, like the one implemented in the geneSetTest R function, and the correlation-corrected version of the test implemented with the CAMERA method.

As I understand it, CAMERA corrects for correlation within a gene set, because a highly correlated query set violates the underlying assumptions and results in a bias toward more significant p-values. So, if I find a significant mean-rank gene set test p-value, but an insignificant CAMERA p-value, I conclude that my query set is correlated, hence the correction factor and reduction in significance.

How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value? Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another? I feel like there's a clear connection between CAMERA's variance inflation factor and a reduction in significance, but I struggle to see how the correction factor can increase the significance over the mean-rank test.

# Answer
> 4 votes
I assume that you are talking about the implementation of these methods in the limma package. Otherwise this answer does not apply.

I think that your questions can be answered with some simulations where we can test with some "genes" with a known relationship:

```
library("limma")
set.seed(123)

# Create some genes and samples
nGenes <- 40
nSamples <- 10
expr <- matrix(rnorm(nGenes*nSamples), ncol = nSamples, nrow = nGenes)
genes <- combn(LETTERS, 2, FUN = paste, collapse = "")
rownames(expr) <- genes[seq_len(nGenes)]
colnames(expr) <- letters[seq_len(nSamples)]

# Create some other genes with slightly higher values
expr2 <- expr+0.5
rownames(expr2) <- genes[seq(41, to = 80)]
expr <- rbind(expr, expr2)

# Some genes are more expressed in the last 5 samples
expr[41:80, 5:10] <- expr[41:80, 5:10] + 2

# Heatmap of the simulated genes and samples
heatmap(expr, scale = "row", main = "Heatmap")

```

```
# Create a design matrix and
design <- matrix(c(rep(1, 10), rep(0, 5), rep(1, 5)),
                 ncol = 2, nrow = 10)
fit <- lmFit(expr, design)
fit2 <- eBayes(fit)
tt <- topTable(fit2, coef = 2, number = Inf, sort.by = "none")

randomGeneSet <- sample(1:80, 5, replace = FALSE)
# Calculate the enrichment with geneSetTest
geneSetTest(1:40, tt$t)
## [1] 1
geneSetTest(41:80, tt$t)
## [1] 9.043842e-13
geneSetTest(randomGeneSet, tt$t)
## [1] 0.9730717

# Explore the interGeneCorrelation and the  variance inflation factor
interGeneCorrelation(expr[1:40, ], design)
## $vif
## [1] 0.2055766
##
## $correlation
## [1] -0.02036983
##
interGeneCorrelation(expr[41:80, ], design)
## $vif
## [1] 12.79611
##
## $correlation
## [1] 0.3024643
##
interGeneCorrelation(expr[randomGeneSet, ], design)
## $vif
## [1] 0.3526011
##
## $correlation
## [1] -0.1618497

# Calculate the enrichment with camera
camera(expr, 1:40, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1     40 -0.02036983      Down 3.0445e-06
camera(expr, 41:80, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction      PValue
## set1     40   0.3024643        Up 0.002402933
camera(expr, randomGeneSet, design, contrast = 2, inter.gene.cor = NA)
##      NGenes Correlation Direction     PValue
## set1      5  -0.1618497      Down 0.04272392

```

In this simulation we created three gene sets were the p-value of both camera and geneSetTest where significant, one where only camera was significant and another one where both were below the 0.05 threshold of significance.

As you can see we can calculate the variance inflation factor for each gene set, which also yields the intercorrelation for each gene set we can summarize it into a table:

```
               camera    VIF    correlation  geneSetTest
GS1 (1:40)     3.0e-06   0.203     -0.020     1.0
GS2 (41:80)    0.024     12.8       0.30      9.0e-13
randomGeneSet  0.043     0.353     -0.16      0.97

```

As GS1 and GS2 are the complete opposite we see that they have quite different inter correlation, which leads to a different VIF. We can see that higher VIF or small doesn't mean that the gene set has lower p-value for camera.

> How can I interpret the opposite case, where I find a significant CAMERA p-value but an insignificant mean-rank gene set test p-value?

I am unsure how to interpret that comparison. It would be better to explore this with real data (for example from ALL library), and contact the authors, either via the support site in Bioconductor or through personal communications.

> Does this imply that my query set is less correlated than expected, or that non-query genes are more highly correlated with one another?

As you might have noticed I used `inter.gene.cor = NA` in order to estimate the correlations and avoid guessing the correlation between the genes but still the correction of intercorrelations has lead to higher p-values. So no, a wrong estimation of the correlation of the query set doesn't explain why the p-value is lower.

Note: I only used the unranked method of the limma 3.34.1, I didn't use the preranked camera (`cameraPR`) because it doesn't have the same results as in `camera` and it doesn't allow to estimate the intercorrelation.

---
Tags: gene, statistics, differential-expression, pathway, gsea
---


title: Coordinates to Genome R Studio Script?

text:
# Question
Title: Coordinates to Genome R Studio Script?
I have a very, very long list of genome coordinates of areas. I need to gather a sequence in order to compare it to another database of sequences that I have. How do I produce the genome sequences from these coordinates in R?

I am very new to this and need a detailed answer.

# Answer
> 2 votes
This doesn't use R, but a simple and quite convenient method is to use `bedtools getfasta` with your fasta file. You'll need to make a BED file out of the RepeatMasker output:

```
awk 'BEGIN{OFS="\t"}{name=sprintf("%s_%s", $10, $11); strand="+"; if($9 == "C"){strand="-"}; print $5,$6-1,$7,name,".", strand}' repeatmasker.output > rmsk.bed

```

Then, given that and your `genome.fa` fasta:

```
bedtools getfasta -name -s -fi genome.fa -bed rmsk.bed -fo output.fa

```

This will produce a multiline fasta file names like `AluSx_SINE/Alu` (due to the `-name` option) where the sequences are reverse complemented if appropriate.

You can install bedtools either via conda (`conda install -c bioconda bedtools`) or github.

---
Tags: r, genome, sequence-alignment, ngs
---


# Question
Title: Canu assembly not making a single consensus?
I've downloaded reads from this BioProject. Using canu with default parameters (no correction), I've got 4 contigs, none of which really look like the reference plasmid here.

The command I used was:

```
canu -p ip40a -d ip40a_assembly -useGrid=false -maxThreads=6 -genomeSize=175k -nanopore-raw reads.fastq

```

Granted, the original publication for this sequence used minimap/miniasm but my questions are:

A. Why do I get multiple contigs? Doesn't canu generate a single assembled sequence?

B. This is a plasmid but canu says `suggestedCircular=false`. What's going on here?

# Answer
> 3 votes
Looking that the link you provided only 30% of the reads appear to be e. coli. I would suggest filtering the reads for those that align to e. coli and see where that gets you. Did the original publication describe the method they used for data cleaning?

A:

```
https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?run=SRR3385327

Unidentified reads: 68.82%
Identified reads: 31.18%
cellular organisms: 30%
Bacteria: 30%
Proteobacteria: 28.09%
Gammaproteobacteria: 19.7%
Enterobacterales: 13.3%
Enterobacteriaceae: 12.54%
Vibrionales: 0.01%
Viruses: 1.17%

```

B:

Without a good assembly canu won't be able to determine if the assembled genome is circular. I would ignore this until you have a good example.

# Answer
> 1 votes
The *E. coli* genome is closer to 4.6 Mb in size; if there are a large quantity of host reads in there, then they'll take over the reads and there won't be enough for your target plasmid. By default, it will only try to assemble from a maximum of ~40X coverage of the specified genome size (higher values do not add any more information for assembly).

It's possible that miniasm will take all reads into account, rather than just a portion of reads. You may need to filter out the *E. coli* reads for this assembly to work properly.

For what it's worth, Canu's default parameters *do* include a correction step.

---
Tags: nanopore, fastq, canu
---


# Question
Title: What would be a reasonable healthy data set/control group to compare with a cell line from TCGA?
Assume I want to run a simple machine learning of healthy/non healthy classification, how can I use TCGA where it has only data of people with cancer. What do usually people do for example for study of gene expressions?

# Answer
In the TCGA datasets there is a variable (type) containing the information if the sample was from the tumor or from an adjacent region, which is usually considered as healthy.

As this samples come from the same patient and the same tissues they are used in the analysis for the comparisons and classifications.
> 2 votes

---
Tags: database, cancer
---


# Question
Title: Interpreting 0x200 flag in bwa-mem alignments
I am looking at the `bwamem.h` code in http://github.com/lh3/bwa and found that BWA-MEM will give flag 0x200 to what it calls `MEM_F_SOFTCLIP`:

```
#define MEM_F_PE        0x2
#define MEM_F_NOPAIRING 0x4
#define MEM_F_ALL       0x8
#define MEM_F_NO_MULTI  0x10
#define MEM_F_NO_RESCUE 0x20
#define MEM_F_REF_HDR   0x100
#define MEM_F_SOFTCLIP  0x200
#define MEM_F_SMARTPE   0x400
#define MEM_F_PRIMARY5  0x800
#define MEM_F_KEEP_SUPP_MAPQ 0x1000

```

But looking at the explanation in this website: https://broadinstitute.github.io/picard/explain-flags.html

It labels `0x200` as read fails platform/vendor quality checks.

If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC? Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?

The SAM/BAM specification labels `0x200` the same way as the Picard flag explaining website:

The question remains as to why BWA-MEM uses the `0x200` flag for softclipping notation rather than Fail QC.

# Answer
bwa mem is using exactly the same meaning of the `0x200` flag as every other program, including picard. Don't blindly assume that that entry in the header file relates to the output flags in the SAM output. If you search the code, you'll see that bwa is using this internally as a flag for "this alignment is soft-clipped", which it handles differently for primary and secondary alignments (hard-clipping the latter).

For what it's worth, `bwa mem` (at least the version I just tried) doesn't set the QC failed flag. This behavour is actually pretty common among aligners, since Illumina has changed how they represent QC failure at least once.
> 3 votes

# Answer
The connection between soft clipping and a read level flag is interesting. Soft clipping occurs at a per base level and is denoted in the CIGAR string. The failQC flag is applied to a read and usually comes directly from the sequencer. It is likely that BWA just named their internal var weird.

> Does the failQC tag come from the SAM/BAM specification or is this something IGV came up with?

The failQC flag is part of the Sam Spec.

> If we have reason to believe that soft clipping can happen for some of the alignments, should we just ignore the fact that IGV calls this fail QC?

If you expect soft clipping I would suggest trimming that sequence from the reads prior to alignment. Also you can look at read in the input reads and see if the failQC flag is set prior to alignment to understand where it is getting set. IGV is just displaying what is in the BAM file not making any decisions about those reads.
> 0 votes

---
Tags: bwa, sam
---


title: About CAFA experiment

text:
# Question
Title: About CAFA experiment
I have a question about CAFA1 and CAFA2 experiment.

CAFA is an experiment that compare computational methods predicting protein functions. In the paper discussing CAFA2, there is an image comparing 2 methods (CAFA1 versus CAFA2):

How do I interpret the numbers and colors (dark and light)? What do they mean?

# Answer
**Note** that in the article published that figure is different from the preprint!
See the figure:

As they aren't substantially different in display other than showing the two set of predictions that are comparable between CAFA1 and CAFA2 I explain them below:

The colors are the margins of $F\_{max}$ they refer to the equation found in the section "Comparison between CAFA1 and CAFA2 methods" under methods:

$$\delta(m\_2, m\_1) = \frac{1}{B}\sum\_{b=1}^{B}F\_{max}^{b}(m\_2)-\frac{1}{B}\sum\_{b=1}^{B}F\_{max}^{b}(m\_1)$$

Where B is the 10000 bootstraps they did and $m\_1$ and $m\_2$ are the CAFA1 and CAFA2 methods.

Where $F\_{max}$ is explained in the article (and in the preprint) as:

> A perfect predictor would be characterized with $F\_{max} = 1 $

So, the higher the (positive) difference is, the better precision and recall (See wikipedia) are for those methods. Note that the scale is from -0.1 to 0.1 not from -1 to 1 so the differences are quite small (One could argue that the log10 would be better to see those small differences)

One can appreciate that in the bar plots at the right, where the differences between each dataset are small in the B section, so the color is white, and in the C panels they have brighter color.

The numbers in the boxes are the percentage that a CAFA2 method did better than a method in CAFA1 for another dataset in the "intersection of the two GO snapshots" as noted in the methods section. (But it is also not clear to me)
> 5 votes

---
Tags: figure-explanation
---


# Question
Title: Can we use non base-called fast5 files in poretools?
I run the MinION MinKNOW without the live base-calling option. We know there is Metrichor and Albacore to perform base-calling after this process. However, I have not done any base-calling yet.

My question is: Is it possible to use directly the `fast5` files with `poretools` to extract a fastq file without doing any previous base-calling process? I tried it and I get an empty fastq file.

My `reads` directory only contains the `fast5` folder and I run:

`poretools fastq fast5/ > output.fastq`

Any ideas why do I get an empty file? What is the difference between base-called `fast5` and non-base-called `fast5` files?

Now I am trying to do base-calling with Albacore to see if I get a fastq file.

# Answer
> 4 votes
No, poretools does not do basecalling. The `poretools fastq` command can be used to extract the FASTQ information from the basecalled FAST5 file (via MinKNOW live-basecalling or albacore). Alternatively, both of these basecallers can export a FASTQ file directly if desired.

The difference between a basecalled FAST5 and a nonbasecalled FAST5 is the presence of the FASTQ within the file. This can be directly viewed using an HDF5 viewer, like h5dump from the commandline, or HDFView from the GUI.

# Answer
> 1 votes
Albacore can be used to create a fastq, directly from your fast5 files. As such there is no more reason to use poretools, expect perhaps for getting some metadata. With regards to metadata I would like to add that the sequencing\_summary.txt generated by albacore also contains a lot of metadata.

---
Tags: minion, albacore
---


title: What do 0s in stageR mean?

text:
# Question
Title: What do 0s in stageR mean?
I'm running stageR as part of a pipeline based on DRIMSeq and wanted to know what it means when the transcript level padj is equal to 0, particularly when the gene is non-significant. I'm pretty sure that in some cases if the padj is 0 it means that the chance of it being different due to chance is 0. There are a lot more 0s than I would expect in sites where the gene padj is non-significant. One example, I have a gene padj that's \>0.9 but both of the sites I'm looking at have a padj of 0.

I have a couple of suspicions about why this is: either because it defaults to this if the comparative difference is or it uses 0 as a substitute to NA. Am I thinking along the right lines? Are the 0s real or are they likely to be a quirk of the package?

# Answer
In stageR, we automatically return stage-wise adjusted p-values of zero for transcripts coming from genes that were found to be significant and only have 1 transcript, when testing for differential transcript expression (DTE). Similarly, when testing for differential transcript usage (DTU), we return p-values of zero for transcripts derived from a significant gene consisting of no more than two transcripts (based on the tx2gene object you provide).

Transcript-level p-values for non-significant genes should not occur. We provide NA's for genes with only one transcript for DTU, and by default you should get NA p-values for non-significant genes, since these genes do not proceed further than the screening stage.

Could you provide some more details on the analysis? Which multiple testing correction have you selected?

Would it be possible to provide a reproducible example for your analysis?
> 4 votes

---
Tags: r, rna-seq
---


title: RcppArmadillo installation error

text:
# Question
Title: RcppArmadillo installation error
When I am trying to install the RcppArmadillo Bioconductor library in R, it is giving me the following error:

```
/usr/bin/ld: cannot find -lgfortran
collect2: error: ld returned 1 exit status
make: *** [RcppArmadillo.so] Error 1
ERROR: compilation failed for package ‘RcppArmadillo’
* removing ‘/home/nikita/R/x86_64-pc-linux-gnu-library/3.2/RcppArmadillo’

```

I am using Ubuntu 12.04, R 3.2.5 and `g++-5.4`. Any help would be greatly appreciated.

# Answer
> 1 votes
I found a similar question on stackoverflow and the solution was to create a symlink as described by `kevin` there:

```
sudo ln -s /usr/lib/x86_64-linux-gnu/libgfortran.so.3 /usr/lib/libgfortran.so

```

# Answer
> 2 votes
On Ubuntu 14.04 there is a direct route to installing RcppArmadillo. You could perhaps try the same syntax.

```
sudo apt-get update
sudo apt-get install r-cran-rcpparmadillo

```

---
Tags: r, software-installation
---


# Question
Title: sort a fasta file containing the Oxford Nanopore Technologies (ONT) header by sequencing start_time ascending
I have basecalled ONT reads and converted them to multifasta. The multifasta contains the original ONT headers in this format:

`3ebb56cd-3671-4abd-b1ac-0c759bt068d0 runid=eb6214851489c8e00eb0dcbd00d737f7dddxxxx read=1078 ch=422 start_time=2017-11-21T20:35:34Z`

By checking the sequencing output as `fasta`, the sequences are not ordered by start\_time.

```
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC

```

I learned how to extract header IDs and `time_start` from ONT fastq files by using this answer.

I'd like to have a similar approach but with fasta.

Yet, I would like to create a new fasta file where the sequences are written in an ascending order as they were sequenced (`start_time`) as follows:

```
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTTTCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATG

```

# Answer
Edit: After doing some reading, if it can be assumed that all the timestamps are ISO8601 and of the same representation of ISO861 (there are a few), then they can be sorted lexicographically, directly, without converting to epoch-time integers. So I have edited this answer to add an answer without conversion steps:

---

Here's one approach that uses the timestamps directly:

(1) Extract the start times from the headers and sort them to another file, in ascending order:

```
$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa | sort > starts_iso8601.sorted.txt

```

(2) Use that file of ordered ISO8601-formatted starts to write out FASTA records in your sort order:

```
$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

```

Quite a bit cleaner and faster.

You might separate the timestamp extraction and sort steps into two commands and inspect the intermediate result, just to ensure that you have valid strings. Or use the approach below to convert the timestamps, which will catch bad input.

---

Here's a second way that may add unnecessary work. As an advantage, however, it adds some error checking by converting the ISO8601-formatted timestamp.

If the string you pull out of the header is not a valid timestamp, the first conversion step (ISO8601-to-epoch) will fail and you will then be able to catch invalid input pretty quickly, which can add confidence that your output is correct.

(1) Extract the start times from the headers:

```
$ awk '(/^>/){ match($0, /start_time=(.*)/, a); print a[1]; }' my.fa > starts_iso8601.txt

```

(2) Convert the ISO8601-formatted times to time-since-epoch times:

```
$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"$line" +%s; done < starts_iso8601.txt > starts_epoch.txt

```

(3) Sort the epoch times — which are just integers — to another file, in ascending order:

```
$ sort -n starts_epoch.txt > starts_epoch.sorted.txt

```

(4) Convert sorted epoch starts back to ISO8601 format:

```
$ while IFS='' read -r line || [[ -n "$line" ]]; do date -d"@$line" -u +"%Y-%m-%dT%H:%M:%SZ"; done < starts_epoch.sorted.txt > starts_iso8601.sorted.txt

```

(5) Use that file of ordered starts to write out FASTA records in your sort order:

```
$ while IFS='' read -r line || [[ -n "$line" ]]; do awk -v line=$line -v RS=">" '{ m = match($0,line); if (m) { printf(">%s", $0); } }' my.fa ; done < starts_iso8601.sorted.txt > my_ordered.fa

```

Note: Use `gawk`, `gdate`, `gsort` in place of `awk`, `date` and `sort` if you're using OS X, Homebrew, and GNU core utilities.

I'm showing the steps separately to explain the general procedure.

Many or all of these steps can be glued together in a Unix pipeline to save the time, memory and disk space costs of writing intermediate files.

Also, I'm not doing much error checking here. You might look at cases where headers do not contain `start_time`, or times are not all in ISO8601, or other cases not anticipated here.
> 3 votes

# Answer
I've just made a small tweak to my fastx-sort.pl script to get this working. It works on both fasta and fastq sequences, and should even work on a random mixture of the two.

This script can now do a dictionary sort by the bits of an ID string that match an arbitrary search pattern, as well as the previous ID and length-based sorting:

```
fastx-sort.pl -p '<perl regex pattern>' <fastq file> > sorted.fastq

```

Here's an example with your fasta input, sorting on the nanopore time string:

```
$ cat seq.fa | grep start_time
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' seq.fa | grep start_time
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z

```

Here's a worked example with one of my fastq files for the nanopore time string:

```
$ cat called_all.fastq | grep start_time | head -n 3
@3a151247-6720-4071-945e-7a92c1856e45 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9896 ch=246 start_time=2017-11-27T00:58:42Z
@52e439d6-f14a-4af2-adba-132d240fa967 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8634 ch=511 start_time=2017-11-27T00:58:41Z
@82e3dd50-c395-4ade-ab6b-0eed530d3c3f runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=8312 ch=299 start_time=2017-11-27T00:58:39Z
$ ~/scripts/fastx-sort.pl -p 'start_time=(.*?Z)' called_all.fastq | grep start_time | head -n 10
@195ad363-1d5a-40bf-abed-820c8acce8fc runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=9 ch=144 start_time=2017-11-26T22:42:32Z
@5df2fc6e-53c8-489c-8002-c97f151104bd runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=11 ch=129 start_time=2017-11-26T22:42:37Z
@23940fe6-0eef-418f-9156-d0a26df9bd80 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=22 ch=330 start_time=2017-11-26T22:42:39Z

```

This script now does a stable sort, so you can do things like sorting by channel, then by read number:

```
~/scripts/fastx-sort.pl -n -p 'read=([0-9]+)' called_all.fastq | ~/scripts/fastx-sort.pl -n -p 'ch=([0-9]+)' | grep 'ch='
...
@876327ec-ccb0-44ad-a9fb-e008744d6c15 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=5974 ch=66 start_time=2017-11-27T00:04:52Z
@e3ec5d92-ee8a-43d6-82ed-e359e0539b75 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6117 ch=66 start_time=2017-11-27T00:06:55Z
@410fe169-c8d1-4686-b92f-cbeb0602b2f0 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6280 ch=66 start_time=2017-11-27T00:09:26Z
@d27c2827-1746-4827-9fec-3075e08d3d64 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=6375 ch=66 start_time=2017-11-27T00:10:44Z
@43079baa-8591-4498-bbfd-5f5ad677c011 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=32 ch=69 start_time=2017-11-26T22:43:32Z
@2bd9b317-438d-4e45-958b-7f77df1dd51a runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=170 ch=69 start_time=2017-11-26T22:47:32Z
@e3599944-849a-492c-ad6d-8cc319b7c5d7 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=207 ch=69 start_time=2017-11-26T22:49:11Z
@8b9d80d0-b28d-46bb-a4e9-6bba95a79ff4 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=320 ch=71 start_time=2017-11-26T22:50:33Z
@b964f164-c36c-4f8a-8a1a-1ad196c996a3 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=780 ch=71 start_time=2017-11-26T22:59:14Z
@92497c10-a46a-4a31-bc18-cfa334491ca9 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=937 ch=71 start_time=2017-11-26T23:02:32Z
@2a715984-6e14-4f31-bd47-d36ee5cbc3e8 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2271 ch=71 start_time=2017-11-26T23:24:11Z
@ff2a59ae-4909-4f3e-bea2-b695d521bf5e runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2441 ch=71 start_time=2017-11-26T23:27:01Z
@f106792d-31a4-44f2-a4fe-83f35a42a573 runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2559 ch=71 start_time=2017-11-26T23:28:54Z
@8de1d1ca-c6f7-43f5-bbd6-3c79187c1bbe runid=256e8ca97106d2fef2538e4c078d1780a9ddb95a read=2930 ch=71 start_time=2017-11-26T23:34:33Z
...

```

Additional notes:

* In order to do the sort, the script currently stores the entirety of the input in memory.
> 2 votes

# Answer
My go-to tools for any type of fasta file manipulation are these two simple awk scripts:

* FastaToTbl will convert fasta to "tbl" formato: the ID, a tab and then the sequence, all on the same line.

```
#!/usr/bin/awk -f
{
        if (substr($1,1,1)==">")
      if (NR>1)
                     printf "\n%s\t", substr($0,2,length($0)-1)
      else
       printf "%s\t", substr($0,2,length($0)-1)
            else
                    printf "%s", $0
}END{printf "\n"}

```
* TblToFasta will convert tbl back to fasta:

```
#!/usr/bin/awk -f
{
  sequence=$NF

  ls = length(sequence)
  is = 1
  fld  = 1
  while (fld < NF)
  {
     if (fld == 1){printf ">"}
     printf "%s " , $fld

     if (fld == NF-1)
      {
        printf "\n"
      }
      fld = fld+1
  }
  while (is <= ls)
  {
    printf "%s\n", substr(sequence,is,60)
    is=is+60
  }
}

```

Armed with those and a bit of Perl, you can do:

```
$ FastaToTbl file.fa |
    perl -lne '/start_time=(\S+)/;
               $d=`date -d $1 +%s`;
               $k{$d}=$_;
               END{
                    print "$k{$_}" for sort(keys(%k))
               }' | TblToFasta
>fa9f83ee-ec58-419b-aff7-e40604cb7a34 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=183 start_time=2017-11-21T20:35:09Z
TGGTATGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTATGTGGTAT
GTGGTATGTGGTATGTGGTATGTGGTATG
>371a075a-6701-4187-affc-7aa267c41850 runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=315 start_time=2017-11-21T20:35:57Z
CTCGTTTCAGGTATGCTTCGTTTCAGTTGCCCTCAGGAGCAGATTGTTTATGGGTTACTT
TCCTGTCGCTCTATCTTCAC
>f546fab3-39b8-4efb-9bc7-cf6ab19cf6dd runid=eb6214851489c8e00eb0dcbd00d73214851489c8 read=100 ch=363 start_time=2017-11-21T20:36:34Z
TTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGGTATGTTGG
TATGTTGGTATGTTGGTATGTTGGTATG

```

### Explanation

The perl script will take the tbl format sequences (so everything on one line), and identify the `start_time`: `/start_time=(\S+)/`. The parentheses capture the matched pattern so the start time will now be `$1`. Then, we run the system `date` command giving it the start time as the date and getting the result in seconds since the epoch (`%s`). note that this assumes GNU `date` (the default on Linux systems) which can handle the `-d` option to specify a date. The next step is using `$d`, the date in seconds, as a key for a hash whose values are the input lines: `$k{$d}=$_`. Then, after the whole file has been processed, we print each entry of the hash in sorted order (`print "$k{$_}" for sort(keys(%k))`). Finally, TblToFasta changes back to fasta format.
> 0 votes

---
Tags: fasta, nanopore, ngs
---


# Question
Title: How to trim adapter sequences from GSE65360 in order to map the reads?
I am trying to map single-cell chromatin accessibility data from doi:10.1038/nature14590, obtained using scATAC-seq, to the reference genome. Example paired-end reads are here.

The reads contain various adapter sequences, which the authors "trimmed from FASTQs using custom python scripts" \[ref\]. I am having problems reproducing this trimming. They used Nextera Transposase Adapters

> TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG

and custom barcode adapters which are listed \[here\]. An excerpt from this list:

> Custom Barcodes Adapter 1 (index i5): AATGATACGGCGACCACCGAGATCTACAC**TAGATCGC**TCGTCGGCAGCGTCAGATGTGTAT AATGATACGGCGACCACCGAGATCTACAC**TAACCAAG**TCGTCGGCAGCGTCAGATGTGTAT AATGATACGGCGACCACCGAGATCTACAC**AGGTTGCC**TCGTCGGCAGCGTCAGATGTGTAT Custom Barcodes Adapter 2 (index i7): CAAGCAGAAGACGGCATACGAGAT**ACCCAGCA**GTCTCGTGGGCTCGGAGATGTG CAAGCAGAAGACGGCATACGAGAT**CCCAACCT**GTCTCGTGGGCTCGGAGATGTG CAAGCAGAAGACGGCATACGAGAT**CACCACAC**GTCTCGTGGGCTCGGAGATGTG

Analyzing the FASTA-files with FastQC showed that the reverse complement of those sequences was over-represented in the data and also confirmed Nextera Transposase Adapters.

I tried to trim the sequences using `cutadapt`:

`cutadapt --minimum-length=20 -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATA -A ATACACATCTGACGCTGCCGACGA -A CTGTCTCTTATA -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz SRR1780164_1.fastq.gz SRR1780164_2.fastq.gz`

with

* `CACATCTCCGAGCCCACGAGAC` being the reverse complement of the 2nd fix part of the i7 adapters
* `ATACACATCTGACGCTGCCGACGA` being the reverse complement of the 2nd fix part of the i5 adapters
* `CTGTCTCTTATA` for the Nextera Transposase Adapters

FastQC does not indicate overrepresented sequences in the resulting FASTA-files, however still indicates Nextera Transposase Adapters in the \_1 files. On the other hand, after TrimGalore removes `CTGTCTCTTATA`, Nextra Adapters are not indicated anymore, for some reason.

Furthermore, mapping the reads using `bowtie2`

`bowtie2 -x GRCh38_primary_assembly_index -1 SRR1780164_1_trimmed.fastq.gz -2 SRR1780164_2_trimmed.fastq.gz > SRR1780164.sam`

results in a low alignment rate

```
1352 (2.79%) aligned concordantly exactly 1 time
3420 (7.07%) aligned concordantly >1 times
11.25% overall alignment rate

```

which gives me the impression I did something wrong during adapter trimming.

Any help would be much appreciated.

--

Edit: As I am not sure about the sequences to use for trimming, I also tried `-a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA`, as found in a suggestion for Nextera Adapters, however it did not remove the sequence over-representation in the -1 file. `-a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA` solves adapter-issues according to FastQC, but also does not give good mapping results.

# Answer
First the Nextera adapters and the custom barcode adapters overlap each other

```
Nextera 1                                     TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Barcode 1   GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT

Nextera 2                                    GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Barcode 2    CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG

```

Looking at some read 1s and read 2s you were correct that the adapters appear at the 3' end of the read reverse complimented (The `|>` points to the adapter in the sequence). You can see that `CTGTCTCTTATACACA` would be the ideal sequence to use for the 3' end of the reads. Since the Nextera adapter and the barcode overlap you only need to use the one that is adjacent to the insert.

```
SRR1780164_2.fastq: GCCCGAGGGGACGCTGCCCCACGTGCTCCAGG |> CTGTCTCTTATACACATCTGACGCTGCCGACGATGCACGAAGTG
SRR1780164_1.fastq: GGGCTGGACAGCGGCGGGGACGCGGACCTGCAG |> CTGTCTCTTATACACATCTCCGAGCCCACGAGACACAAAGTGA

```

Of the 64k reads there are ~11k with the Nextera adapter at the 3' end based on a rough sequence `grep`

On the 5' end of the read there is some Nextera adapter present but it is not reverse complimented and there are only a handful of reads with the adapter present ~250 total between read 1 and read 2.

```
SRR1780164_2.fastq: TGGGCTCGGAGATGTGTATAAGAGACAG <| CTCCCCCGCCTCCAGCATCCGGGCGAGGTAGTGCATCGACGCGTCCAC

```

I used this command:

```
cutadapt --minimum-length=20 -A CTGTCTCTTATACACA -a CTGTCTCTTATACACA -G GATGTGTATAAGAGACAG
-g GATGTGTATAAGAGACAG -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz
SRR1780164_1.fastq SRR1780164_2.fastq

```

Which gives a trimming result of

```
   === Summary ===

Total read pairs processed:             64,832
  Read 1 with adapter:                  16,721 (25.8%)
  Read 2 with adapter:                  16,996 (26.2%)
Pairs that were too short:                  46 (0.1%)
Pairs written (passing filters):        64,786 (99.9%)

```

I think this best represents the library construct and the trimming of the adapter sequences
> 4 votes

---
Tags: read-mapping, single-cell, trimming, adapter
---


title: Finding data to validate results

text:
# Question
Title: Finding data to validate results
My group has a complex datset and I would like to validate the results and check the methods I use in a similar dataset (better if it has been already studied).

Characteristics of my dataset:

* It is related to the inflammatory bowel disease (IBD)

More specifically to Crohn's disease.
* It has RNA-seq

I have explored the GEO and there are few projects (10) with RNA-seq on the tissue of my disease, and I don't know how to find the 16S-seq of those projects if it exists (see point below).
* It has 16S-seq (and not only from the stools)

The ENA seems to have 16S sequencing but most of them are from stools not from biopsies and still, I don't know how to find the RNA-seq of those projects if it exists (see point above)
* It has both RNA-seq and 16S-seq from the *same* individuals at the *same* time point.

The NIH Human Microbiome Project for my disease has only one individual with that both 16S-seq and RNa-seq of the same patient at the same time.

How can I find a dataset similar to mine?

# Answer
> 2 votes
The snarky answer is that it's better to find a matching dataset at the experimental design phase (i.e. *before* these studies are carried out), so that the bioinformatics component on the other side of things ends up being a lot easier and more robust. Although it can be quite time consuming to do a database search at the start, it's even harder to find matching projects if the legwork hasn't been done during the design phase. Researchers should be thinking about how hypotheses will be validated *before* the experiments are carried out.

Every additional characteristic that is chosen for comparison adds complications with regards to matching data, and reduces the applicability of results when there is a mismatch. You would have to be somewhat lucky to be able to find a IBD-Crohn's dataset with paired RNA-seq data and 16s-seq (presumably short-read; is 454 good enough?) from a specific non-stool tissue.

Apart from the places you've already looked, there's dbGaP and SRA. Given the number of characteristics that have been mentioned, I expect you'll have to make some compromises when comparing other datasets.

---
Tags: human
---


# Question
Title: Link Prediction in Bipartite Networks between biological pathways and drugs
Currently I'm building a big matrix (using microarray, mass spectrometry, RNAseq data) that consist by pathways (rows) and treatment/drugs (columns).

The values of that matrix are scores that describe the effect/relation of that specific treatment/drug on each pathway.

Now, some of the drugs is known to target specific disease (target disease) and thus is expected to have a high score with that disease's pathways. **(A)**

The rest of the drugs is known to target other diseases.**(B)**

My thought now is that, if I label the known drugs as "1" (or whatever), could I use a prediction model to find out new relationships (links) between drugs from the **B** class and the significant disease pathway from **A** ?

The whole thing, we can represent it as a bipartite network shown below.

For example, here the denoted as green, drugs are the known drugs belong to the class **A** that mentioned before and they are connected highly with some pathways of the target disease.

From the other hand blue and red ones are drugs taken from other cases (class **B**) that are used to target other diseases, but as my plot shows, blue seems to has an "affinity" to some of the pathways of the target disease and thus should be an alternative candidate for that disease. The red from the other side, doesn't have any "affinity" to any significant pathway and thus we can not consider it as possible candidate.

To find a solution to that kind of problem, I looked for algorithms and methods developed for social networks (e.g recommendation etc.) and mostly I read about random walks with return.

But before to proceed, I would like to know if someone else tried such an approach here and having more experience.

# Answer
> 2 votes
I am not sure if I understood correctly how you classify the drugs. But what you attempt to do is similar to what they do here but in that article they use the drugs targeted to a molecule/pathway to find combination of drugs that are better for the disease (not an alternative drug), see the first image:

In that article they mention the pathway cross-talk inhibition (PCI), you might want to calculate it for your pathways to assess the relative reduction of the pathway efficiency.

---
Tags: networks, pathway, machine-learning, modelling
---


# Question
Title: How can I apply proportional (p) distances (Nucleotide) using bioPython
I am working on phylogenetic tree generation. I used BioPython and ClustalW 2.x.x for this purpose. I have generate the tree using BioPython but when I try to generate tree using "MegaSoftware GUI" tool my tree does not matched with the output of MegaSoftware.

When I use "No. of difference" model then results are matched with biopython tree Like below image shown.

But whenever I choose "p-distance" model my results are not matched.

**Output MegaSoftware using p-distance**

**Output BioPython**

Here is my code for generating phylogenetic tree:

```
def conversion_fasta_dnd(clustalw, filepath, type, clustering, output, pwmatrix):
    """
        ClustalW tool to conversion the file format.
        :param clustalw: path of the executable programm.
        :param filepath: path of `fasta` extension file.
        :param type: dna, protein
        :param clustering: upgma, nj
        :param output: clustal, phylip, fasta
        :param pwmatrix: blosum, pam, id
        :return: dnd, aln (files)
        """
    # esblish the full path to the clustalW program
    clustalw_exe = clustalw

    # check to make sure the program is there
    assert os.path.isfile(clustalw_exe), "Clustal W executable missing"

    # create the appropriate command line
    clustalw_cline = ClustalwCommandline(clustalw_exe,
                                         infile=filepath,
                                         type=type,
                                         clustering=clustering,
                                         output=output,
                                         pwmatrix=pwmatrix)

    # print the command line
    print(clustalw_cline)

    # execute the command
    stdout, stderr = clustalw_cline()

# This function draw the tree on image file
def simple_tree(filepath):
    """
    Generate the simple tree image.
    :param filepath: path of dnd file format.
    """
    def get_label(leaf):
        return leaf.name

    file = filepath
    tree = Phylo.read(file, 'newick')
    Phylo.draw(tree, label_func=get_label, do_show=False, branch_labels=lambda c: c.branch_length)

    # axis management
    pylab.axis('on')
    pylab.savefig(OUTPUT_IMAGE, format='png', bbox_inches='tight', dpi=100, bbox_extra_artists=[])

```

In above code I have used two things first BioPhython lib and the other is Clustal W for align the sequences.

**Could you please clarify what p-distances are?**
I only have these reference for p-distance:

I need to add p-distance model in my code, someone please help me out How can I do that.

# Answer
Based on the sites you mention, it seems that p-distances are pairwise distances between aligned sequences obtained by counting the number of differences and dividing by the length of the alignment. This is the proportion of differences between two sequences.

I'm not sure exactly what you mean by "add p-distance model in my code", but here is how you could compute those distances from an existing alignment (here I use as an example an alignment I have in fasta format):

```
from Bio import AlignIO

def p_dist(seq1, seq2):
    ali_len = len(seq1)
    assert ali_len == len(seq2)
    return sum(0 if nuc1 == nuc2 else 1 for (nuc1, nuc2) in zip(seq1, seq2)) / ali_len

ali = AlignIO.read("most_targeted_histones/transcript_alignments/H4_top4_loss_spliced_NT.fa", "fasta")
nb_seq = len(ali)
distances = {}
for i in range(nb_seq - 1):
    rec1 = ali[i]
    j = i + 1
    while j < nb_seq:
        rec2 = ali[j]
        distances[(rec1.id, rec2.id)] = p_dist(rec1.seq, rec2.seq)
        j += 1

print(*distances.items(), sep="\n")

```

Output:

```
(('ZK131.4_spliced', 'ZK131.8_spliced'), 0.0)
(('ZK131.4_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.4_spliced', 'ZK131.1_spliced'), 0.0)
(('ZK131.8_spliced', 'T10C6.14_spliced'), 0.10897435897435898)
(('ZK131.8_spliced', 'ZK131.1_spliced'), 0.0)
(('T10C6.14_spliced', 'ZK131.1_spliced'), 0.10897435897435898)

```

If you want to build a tree from those distances, I'm not sure you can do this with BioPython. It seems however that there are some possibilities with DendroPy.

You would need to write the distances in a suitable format, and to choose the method between UPGMA and neighbour-joining.
> 3 votes

---
Tags: python
---


# Question
Title: At what processing step should library strandedness type be taken into account?
Suppose I have single-end RNA-seq data for which the reads in the fastq file are reversed with respect to the original extracted RNAs.

Suppose I have the following workflow:

1. Map the reads on the reference genome, sort and index the resulting sam output.
2. Use the bam file to make a bigwig file showing reads mapping on the plus and minus strands of the chromosomes.
3. Use the bam file and an annotation file to quantify the gene expression.

Steps 2 and 3 depend on step 1, but are not dependant on one another.

I'm sure that strandedness needs to be taken into account no later than at step 3.

The questions are the following:

* Is it common to have mapping tools with an option to "correct" for strandedness?
* Assuming the bam file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Regarding this second point, I'm split between two attitudes:

1. There should be no correction, the bigwig files are there to represent the "raw" mapping results.
2. The bigwig files are there to give an idea of the origin and abundance of the originally extracted RNAs, so a strandedness correction should be used if necessary.

Which one is the most common attitude?

If the second one is preferred, then wouldn't it be "safer" / "cleaner" to handle strandedness as early as possible? Possibly at the mapping, or even using a pre-processing step?

# Answer
I usually take into account library strandness during the mapping to the reference step (but I never worked with bigwig files myself).

Assuming RNA-seq data from Illumina, you can use Hisat2 for alignment to the reference genome, and you have the option `--rna-strandness` to specify the strand-specificity information (default unstranded). This is particularly important when you wish to have XS strand attribute in your bam file for easier compatibility with downstream software (especially important in my experience if you are planning to use some utility from Cufflinks for your step 3).
> 2 votes

# Answer
1. If your aligner has a strandedness option then go ahead and use it. The general idea here being that you'll preferentially align correctly to genes.
2. Whether to bother with strandedness here depends on your goal. Are you interested in anti-sense transcription? Do you need to use the bigWig track to accurately quantify translational pausing or some other strand-specific phenomenon? If yes, then go ahead and make different bigWig files for each strand. If you just need the bigWig file for convenience in IGV or creating heatmaps (e.g., with deepTools) then you'll generally be fine without accounting for strand. For well-annotated species, the rate of anti-sense transcription is generally low and the rate of genes overlapping isn't absurdly high. Consequently, you can usually get away with using an unstranded bigWig.
3. If you're using featureCounts you'll certainly want to use the appropriate strandedness setting.

It's relatively common for RNAseq aligners/mappers to offer strandedness options. Certainly tophat2 and hisat(2) do, but also things like salmon.

Again, for bigWig files the question becomes whether you really need two files to get the information you need. You're not going to delete the BAM files, so you can always make them if you need them.
> 2 votes

# Answer
Firstly, right at the start if the experiment it's important that your RNA samples are processed with a strand-specific protocol (e.g. Illumina's TruSeq Stranded) in order to produce stranded libraries for sequencing. If the samples haven't been treated as such, then there's nothing you can do to correct it.

Secondly, I'm not sure stranded protocols are compatible with single-end reads. The aligners I'm aware of use the direction of aligned read pairs to assess the strand from which the RNA fragment has come from. How that works with single-end reads, I don't know?

Finally, to answer your specific question, essentially step 1 is where you need to start with strand-specific options and do so at all subsequent steps as appropriate.
> 1 votes

# Answer
Strandedness should be taken into account at the mapping stage.

If you are presented with only a BAM file, ask for the raw reads. If that's not possible, you can recover the reads (as fastq sequences) from the BAM file and/or remap them, but that conversion requires a lot more effort (e.g. most mappers expect the fastq files from each end to be sorted by read name and in the same order).

Regarding your specific questions...

> Is it common to have mapping tools with an option to "correct" for strandedness?

Yes, mapping tools can be provided with strand direction information for mapping. The way that this is done depends on the program, for example bowtie2 has options `--fr`, `--rf`, and `--ff` to take strandedness into account. I haven't done any of my own comparisons to determine whether these options actually make a difference in mapping for Bowtie2. I know that Salmon (which does mapping and quantification at the same time) cares a lot more about strandedness, and will complain/warn if a stranded library type is discovered and not specified.

> Assuming the BAM file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Yes, but it'd be better to remap with \[the correct\] strandedness information included. When I generate bigwig files to demonstrate strandedness (even assuming that strandedness is incorporated into the mapping), I still need to do a fair amount of legwork to convert that into bigwig files. The data needs to first be split into four groups (e.g. with `samtools view -F 0xXX -f 0xXX`):

1. First read, mapped in the reverse orientation
2. Second read, mapped in the reverse orientation
3. First read, mapped in the forward orientation
4. Second read, mapped in the forward orientation

I then combine two of those groups, depending on the library strandedness scheme. For a typical stranded Illumina run, groups 1 and 4 are combined (e.g. with `samtools merge`), as are groups 2 and 3. These combined groups are turned into separate BigWig files, representing the "genome-oriented" and "genome-reversed" read pairs.
> 1 votes

---
Tags: read-mapping, strandedness
---


title: Identification of Microorganism

text:
# Question
Title: Identification of Microorganism
I have the 16S DNA Sequence of a Microorganism. Can I confirm its identity by only running a BLAST or do I have to follow other methods. If so please mention the methods.

# Answer
As always with bioinformatics, any reasonable answer will be dependendant on what you want to do with this sequence, and you haven't provided enough information in that regard.

For most exploratory analysis BLAST is good enough. You don't *have to* follow any existing methods. The benefit of using NCBI's web BLAST server is that it provides you with a very quick way to find a match from the sequence information in GenBank.

BLAST does a local alignment. If you have the full 16s sequence, and it's a novel organism, you *might* get different match results from a program that does global alignment, but then again you might get different results from a different local alignment tool.
> 3 votes

---
Tags: ngs
---


title: Correlation network representation

text:
# Question
Title: Correlation network representation
I have troubles understanding the figure 10 of this paper:

> Similarity network of participating methods for BPO. Similarities are computed as Pearson’s correlation coefficient between methods, with a 0.75 cutoff for illustration purposes. A unique color is assigned to all methods submitted under the same principal investigator. Not evaluated (organizers’) methods are shown in triangles, while benchmark methods (Naïve and BLAST) are shown in squares. The top-ten methods are highlighted with enlarged nodes and circled in red. The edge width indicates the strength of similarity. Nodes are labeled with the name of the methods followed by “-team(model)” if multiple teams/models were submitted

This study calculated the pairwise Pearson correlation between methods on a common set of gene-concept pairs and then visualized these similarities as networks. How is the place of the methods determined in this visualization?

For example, the relationship between Pannzer (1) and Pannzer (2) is greater than the relationship between Pannzer (1) and Pannzer (3), but the Panzer (1) and Panzer (3) are closer to each other than Panzer (1) and Panzer (2). Is this proximity positioning random? or is there a certain rule?

# Answer
When plotting networks where only the weight of an edge is known like this one. The relative position between elements of the networks is irrelevant. The only informative part of the graph is (usually) the size of the edges, where is the thicker means more weight. This isn't a PCA, so the distance between different nodes is not relevant.

For those reasons there are several methods to build these networks, some of them give more emphasis to represent better the degree of the network or other characteristics of the networks.

---

I consider this representations bad for informing about the networks, because the cut off is arbitrary and the distance between nodes is not informative. One could first do a PCA or MDS and then add the edges between the nodes, this way the relative distance between the nodes would be meaningful as well as the edges. Or one could also use a heatmap to represent all the correlations between the nodes.
> 1 votes

# Answer
The description you quote only states that edge thickness has relevance, not edge length. Usually with this sort of image, the program generating the graph (most probably Cytoscape) will simply arrange the nodes to maximize legibility.

While the edge length can, sometimes, be relevant it usually isn't and instead it is the thickness of the edges and the size of the nodes that carry the information. The positioning of the nodes is most likely arranged to make the image as clear as possible.
> 0 votes

---
Tags: figure-explanation
---


# Question
Title: RapMap: reference transcriptome for simulated reads
I am very new to bioinformatics and trying to repeat the benchmark in the RapMap paper with an experimental tool working in a similar but different fashion.

In the paper (taken from their github) about RapMap, an alignment/mapping tool for rna-seq reads to a reference genome, the dataset used for their benchmarks is described like this:

> To test the accuracy of different mapping and alignment tools in a scenario where we know the true origin of each read, we generated data using the Flux Simulator (Griebel et al., 2012). This synthetic dataset was generated for the human transcriptome from an annotation taken from the ENSEMBL (Cunningham et al., 2015) database consisting of 86090 transcripts corresponding to protein-coding genes. The dataset consists of 48 million 76 bp, paired-end reads. \[...\] When benchmarking these methods, reads were aligned directly to the transcriptome, rather than to the genome.

I understand that they generated reads using flux simulator, but I don't understand which transcriptome they mapped them against. I can't find any other information in the paper about which annotation they used either.

I am assuming one would need to generate an entire transcriptome in .fa format from the GTF and genome? Or did they generate a transcriptome from only the simulated reads and then mapped the same reads against it?

GENCODE offers GTFs and FASTAs for gene-encoding transcripts here, might that be the kind of thing they used? Would the reads simulated from these GTFs fit to these transcript FASTA files?

# Answer
> 5 votes
Thanks for your interest in RapMap.

At that time we were using flux simulator for simulating read sequence data. We used the genome and gtf file together as an input to flux.

I dug into the scripts and got hold of the gtf link `ftp://ftp.ensembl.org/pub/release-80/gtf/homo_sapiens/Homo_sapiens.GRCh38.80.gtf.gz`, and the genome file link `ftp://ftp.ensembl.org/pub/release-80/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.toplevel.fa`.

We used `rsem-prepare-reference` to create the reference transcriptome file later in order to run RapMap and other transcriptome based quantification tools. I hope this helps.

PS: I am a coauthor of the RapMap paper.

# Answer
> 1 votes
If they generated reads from the transcriptome then the started with the FASTA version of the ENSEMBL transcriptome and then aligned the reads back to it. For a project like this I don't think that you need the GTF files, since you don't need to know the location of the transcript on the larger genomic reference. If you can't find an ENSEMBL transcriptome then you can use the GTF and the bedtools program to create a FASTA.

---
Tags: sequence-alignment, read-mapping, transcriptome
---


# Question
Title: HMM Profile from Convergence of retrotransposons in oomycetes and plants
I am trying to recreate the experiment from *Convergence of retrotransposons in oomycetes and plants* by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. In the "Methods" section a

> HMM profile was constructed from the amino acid alignment of aRNH sequences from the Ustyantsev et al

I was wondering where I could get this profile. I checked the cited papers but they were also referencing other papers. Could someone maybe clarify where I could get the profile?

# Answer
> 2 votes
The alignment they used to build the HMM is included as supplementary file 8 in the paper you cite:

> Multiple amino acid sequence alignment of RNH genes and domains from diverse taxa constructed and used for the phylogenetic reconstruction in the present study. (TXT 45 kb)

So, you can download that alignment and then use a tool like hmmbuild (part of HMMER) to build the HMM profile:

```
hmmbuild model.hmm 13100_2017_87_MOESM8_ESM.txt

```

The command above will create the file `model.hmm` which contains the HMM built from the alignment in the supplementary file `13100_2017_87_MOESM8_ESM.txt`. I am including the first few lines here (this is *not* the entire file, don't use this as the model—build it yourself. I only include it to illustrate the format):

```
HMMER3/f [3.1b2 | February 2015]
NAME  13100_2017_87_MOESM8_ESM
LENG  121
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Mon Dec 11 13:33:29 2017
NSEQ  165
EFFN  5.697556
CKSUM 2893352190
STATS LOCAL MSV       -9.8975  0.71332
STATS LOCAL VITERBI  -10.6450  0.71332
STATS LOCAL FORWARD   -3.9247  0.71332
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.34473  4.21240  2.85435  2.69203  3.54228  2.79153  3.54812  3.01140  2.72630  2.48370  3.82225  2.94269  3.49156  3.05882  2.84632  2.63343  2.86306  2.72633  4.27012  3.56181
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.13773  5.45955  2.08419  0.61958  0.77255  0.00000        *
      1   2.99488  4.57852  4.44659  3.29307  1.80597  4.16203  3.88788  2.87265  3.75733  2.38594  3.30344  3.33633  4.53284  3.64324  3.69630  2.90350  3.30742  2.83413  4.36991  1.50673      1 y - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00723  5.32905  6.05140  0.61958  0.77255  0.45275  1.01027
      2   2.66602  2.17750  5.01715  4.40788  1.72617  4.29647  4.62737  2.89579  4.21116  2.71668  3.64495  4.43195  4.66188  4.35829  3.99238  3.10363  2.02207  1.89528  5.14692  3.02155      2 f - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.00661  5.41881  6.14115  0.61958  0.77255  0.61842  0.77391

```

---
Tags: hidden-markov-models
---


# Question
Title: What are the differences between GWAS and QTL mapping?
What are the differences between the two methods?

What advantages does one have over the other, and what are their limitations relative to one another?

# Answer
> 2 votes
Most papers that I've looked at distinguish between looking for genotype/trait associations in mapping populations created by experimental crosses (QTL mapping) and looking for genotype/trait associations in populations of unrelated individuals (GWAS). No distinction is made between susceptibility to disease (often referred to as a case-control GWAS) and continuous traits (all of the big studies on genetic associations with height use the term GWAS, for example).

Many of the papers in my first link use both methods, so you can take a look at them to get an idea of the advantages and limitations of each.

---
Tags: genomics, gwas
---


# Question
Title: How to use the hmmsearch on prototypic repetitive sequences Repbase Update database?
I want to use the hmmsearch proposed in *Convergence of retrotransposons in oomycetes and plants* by Kirill Ustyantsev, Alexandr Blinov and Georgy Smyshlyaev. After the help from terdon I managed to create the correct profile I assumed that the command

```
hmmsearch aRNH_profile.hmm prirep.ref > output.out

```

would be the correct command. I downloaded the compressed database `RepBase22.11.fasta.tar.gz` from here in the fasta format. However when unpacking I get the following tree

```
.
├── angrep.ref
├── appendix
│   ├── angapp.ref
│   ├── athapp.ref
│   ├── bctapp.ref
│   ├── cbrapp.ref
│  ....
├── athrep.ref
├── bctrep.ref
├── cbrrep.ref
...

```

which file is the correct fasta file? I assumed `prirep.ref` might be correct file but the output shows

```
[No hits detected that satisfy reporting thresholds]

```

# Answer
> 1 votes
RepBase explain what these files are on their website. It is usually a good idea to check the documentation when downloading files so you know what you're downloading.

Based on the list of files described on the linked page, there is no `prirep.ref` file, so I assume you mean `plnrep.ref` which is "Other Plants". In any case, again according to that page, the files with plant and/or fungal data (which is what I assume you want; you never said) are:

athrep.ref Arabidopsis thaliana oryrep.ref Oryza sativa dcotrep.ref Other Dicot Plants mcotrep.ref Other Monocot Plants grasrep.ref Other Grass Plants plnrep.ref Other Plants fngrep.ref Fungi

There's also `chlrep.ref` for the green alga *Chlamydomonas reinhardtii*, and `diarep.ref` for the diatom *Thalassiosira pseudonana*. I don't know whether you would want to consider either here.

In any case, since these are fasta files, you cna combine them by simply `cat`ing them together:

```
cat athrep.ref oryrep.ref dcotrep.ref mcotrep.ref grasrep.ref \
    plnrep.ref fngrep.ref > plantsAndFungi.ref

```

The, run your search on the `plantsAndFungi.ref` file.

---
Tags: fasta, hidden-markov-models
---


# Question
Title: low-memory high-speed replacement for Picard MarkDuplicates
I am running Picard MarkDuplicates with the following parameters below. On the file described, it takes about 41.6Gb of RAM memory and about 20-25 minutes to compute (only uses 1 core AFAICS).

```
java -Xmx54890m -jar picard.jar MarkDuplicates \
    I=/home/dnanexus/in/sorted_bam/TST73-86-2IC_S16_L00.bml.GRCh38.karyo.bam \
    O=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.deduplicated.bam \
    M=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.duplication_metrics \
    CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT REMOVE_DUPLICATES=true

```

What equivalent free/open source software (samtools?, sambamba?) can I use to obtain the same deduplicated bam file but at a lower cost, considering the instance used here has 8 cores but only one of them is being used, and has about 68Gb of RAM out of which it's using 46Gb odd at peak.

I think previous versions of samtools/sambamba didn't behave the same way as Picard MarkDuplicates with the parameters above, but I would like to know which recent versions would give me the same output.

# Answer
> 2 votes
As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.

---
Tags: sam, samtools, sambamba, picard
---


# Question
Title: AND or OR gene - protein logical relation in KEGG
How do I know if it is a AND or OR relationship between genes associated with the same enzyme in KEGG? E.g., for EC:1.6.1.2 (http://www.genome.jp/dbget-bin/www\_bget?ec:1.6.1.2) does the bacterium *Streptomyces coelicolor* (SCO) need SCO7622 AND SCO7623 or SCO7622 OR SCO7623?

# Answer
> 1 votes
The enzyme EC:1.6.1.2 is an NAD(P) transhydrogenase and consists of two subunits. The two subunits are encoded by two different genes pntA and pntB.

So in terms of computer logic I would say AND, both subunits are necessary for a functional NAD(P) transhydrogenase enzyme (EC:1.6.1.2).

---
Tags: gene, kegg
---


# Question
Title: NCBI - edirect suite to download all genome sequences associated with a query - troubleshooting
Not sure if this is allowed but I cannot think of a better place to ask. I am attempting to download all genomic sequences from refseq associated with the query Peptostreptococcaceae. I have been been using the following command (based on some reading of the manual):

```
esearch -db nuccore -query "Peptostreptococcaceae" | efilter -molecule genomic -source refseq | efetch -format fasta

```

But I keep getting the following error

```
Unknown option: molecule
Unknown option: source
Must supply -query or -days or -mindate and -maxdate arguments on command
line Db value not found in fetch input No do_post output returned from
'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi?db=nuccore&term=Peptostreptococcaceae&retmax=0&usehistory=y&edirect=3.80&tool=edirect&email=lgunn@sonic15' Result of do_post http request is .............
(a lot of text which I assume is .html or .xml related)

```

# Answer
> 1 votes
The issue here was I was using an older version of eDirect. There is no command to check version or updates. Your best bet is if you experience problems uninstall and reinstall.

# Answer
> 1 votes
If you want to find version info you can either :

Ask for another edirect tool version (it will be the same version for edirect) :

```
efetch -version

```

> 6.80

Check the info contained in the edirect.pl script :

```
grep -i 'version = ' edirect.pl

```

> $version = "6.80";

---
Tags: data-retrieval, edirect
---