title: Cromwell: can't find file ID

text:
# Question
Title: Cromwell: can't find file ID
I am using Cromwell to test a CWL (common workflow language) implementation of the GATK somatic CNV (copy number variation) pipeline. I am currently getting the following error (full log at the end):

> Cannot find a tool or workflow with ID file:///var/folders/cm/0j97xqy577z2dlfht90hv\_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl in file /var/folders/cm/0j97xqy577z2dlfht90hv\_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl

The ID's are the same except for the first starts with "file://". I am not very familiar with Cromwell so I'm not sure how to proceed in resolving this error. Any help would be appreciated.

Thank you

Full log:

```

[2020-12-03 15:20:24,82] [info] Running with database db.url = jdbc:hsqldb:mem:28b7f7a3-d261-42ee-b941-e329cdcf46ba;shutdown=false;hsqldb.tx=mvcc
[2020-12-03 15:20:29,97] [info] Running migration RenameWorkflowOptionsInMetadata with a read batch size of 100000 and a write batch size of 100000
[2020-12-03 15:20:29,98] [info] [RenameWorkflowOptionsInMetadata] 100%
[2020-12-03 15:20:30,03] [info] Running with database db.url = jdbc:hsqldb:mem:ae8f4f13-d319-4439-8681-b684e13ca6ac;shutdown=false;hsqldb.tx=mvcc
[2020-12-03 15:20:30,23] [info] Slf4jLogger started
[2020-12-03 15:20:30,25] [info] Pre Processing Workflow...
[2020-12-03 15:20:30,27] [info] Pre Processing Inputs...
[2020-12-03 15:20:30,46] [info] Workflow heartbeat configuration:
{
"cromwellId" : "cromid-d8084bc",
"heartbeatInterval" : "2 minutes",
"ttl" : "10 minutes",
"writeBatchSize" : 10000,
"writeThreshold" : 10000
}
[2020-12-03 15:20:30,48] [info] Metadata summary refreshing every 2 seconds.
[2020-12-03 15:20:30,50] [info] KvWriteActor configured to flush with batch size 200 and process rate 5 seconds.
[2020-12-03 15:20:30,50] [info] WriteMetadataActor configured to flush with batch size 200 and process rate 5 seconds.
[2020-12-03 15:20:30,50] [info] CallCacheWriteActor configured to flush with batch size 100 and process rate 3 seconds.
[2020-12-03 15:20:30,78] [info] JobExecutionTokenDispenser - Distribution rate: 50 per 1 seconds.
[2020-12-03 15:20:30,80] [info] SingleWorkflowRunnerActor: Submitting workflow
[2020-12-03 15:20:30,82] [info] cwl (Unspecified version) workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a submitted
[2020-12-03 15:20:30,82] [info] SingleWorkflowRunnerActor: Workflow submitted [38;5;2m837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a[0m
[2020-12-03 15:20:30,82] [info] 1 new workflows fetched
[2020-12-03 15:20:30,82] [info] WorkflowManagerActor Starting workflow [38;5;2m837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a[0m
[2020-12-03 15:20:30,83] [[38;5;220mwarn[0m] SingleWorkflowRunnerActor: received unexpected message: Done in state RunningSwraData
[2020-12-03 15:20:30,84] [[38;5;220mwarn[0m] Couldn't find a suitable DSN, defaulting to a Noop one.
[2020-12-03 15:20:30,84] [info] Using noop to send events.
[2020-12-03 15:20:30,85] [info] WorkflowManagerActor Successfully started WorkflowActor-837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a
[2020-12-03 15:20:30,85] [info] Retrieved 1 workflows from the WorkflowStoreActor
[2020-12-03 15:20:30,85] [info] WorkflowStoreHeartbeatWriteActor configured to flush with batch size 10000 and process rate 2 minutes.
[2020-12-03 15:20:30,86] [info] MaterializeWorkflowDescriptorActor [[38;5;2m837d7e7d[0m]: Parsing workflow as CWL v1.0
[2020-12-03 15:20:30,91] [info] Pre-Processing /var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl
[2020-12-03 15:20:32,24] [[38;5;1merror[0m] WorkflowManagerActor Workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a failed (during MaterializingWorkflowDescriptorState): cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anon$1: Workflow input processing failed:
Cannot find a tool or workflow with ID file:///var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl in file /var/folders/cm/0j97xqy577z2dlfht90hv_6h0000gq/T/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.temp.8067760554642643818/837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a.cwl
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.cromwell$engine$workflow$lifecycle$materialization$MaterializeWorkflowDescriptorActor$$workflowInitializationFailed(MaterializeWorkflowDescriptorActor.scala:200)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anonfun$2.applyOrElse(MaterializeWorkflowDescriptorActor.scala:170)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor$$anonfun$2.applyOrElse(MaterializeWorkflowDescriptorActor.scala:165)
at scala.runtime.AbstractPartialFunction.apply(AbstractPartialFunction.scala:34)
at akka.actor.FSM.processEvent(FSM.scala:665)
at akka.actor.FSM.processEvent$(FSM.scala:662)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.akka$actor$LoggingFSM$$super$processEvent(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.LoggingFSM.processEvent(FSM.scala:801)
at akka.actor.LoggingFSM.processEvent$(FSM.scala:783)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.processEvent(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.FSM.akka$actor$FSM$$processMsg(FSM.scala:659)
at akka.actor.FSM$$anonfun$receive$1.applyOrElse(FSM.scala:653)
at akka.actor.Actor.aroundReceive(Actor.scala:514)
at akka.actor.Actor.aroundReceive$(Actor.scala:512)
at cromwell.engine.workflow.lifecycle.materialization.MaterializeWorkflowDescriptorActor.aroundReceive(MaterializeWorkflowDescriptorActor.scala:123)
at akka.actor.ActorCell.receiveMessage(ActorCell.scala:527)
at akka.actor.ActorCell.invoke(ActorCell.scala:496)
at akka.dispatch.Mailbox.processMailbox(Mailbox.scala:257)
at akka.dispatch.Mailbox.run(Mailbox.scala:224)
at akka.dispatch.Mailbox.exec(Mailbox.scala:234)
at akka.dispatch.forkjoin.ForkJoinTask.doExec(ForkJoinTask.java:260)
at akka.dispatch.forkjoin.ForkJoinPool$WorkQueue.runTask(ForkJoinPool.java:1339)
at akka.dispatch.forkjoin.ForkJoinPool.runWorker(ForkJoinPool.java:1979)
at akka.dispatch.forkjoin.ForkJoinWorkerThread.run(ForkJoinWorkerThread.java:107)

[2020-12-03 15:20:32,24] [info] WorkflowManagerActor WorkflowActor-837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a is in a terminal state: WorkflowFailedState
[2020-12-03 15:20:37,26] [info] SingleWorkflowRunnerActor workflow finished with status 'Failed'.
[2020-12-03 15:20:40,51] [info] Workflow polling stopped
[2020-12-03 15:20:40,52] [info] Shutting down WorkflowStoreActor - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Shutting down WorkflowLogCopyRouter - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Shutting down JobExecutionTokenDispenser - Timeout = 5 seconds
[2020-12-03 15:20:40,52] [info] Aborting all running workflows.
[2020-12-03 15:20:40,52] [info] JobExecutionTokenDispenser stopped
[2020-12-03 15:20:40,52] [info] WorkflowStoreActor stopped
[2020-12-03 15:20:40,53] [info] WorkflowLogCopyRouter stopped
[2020-12-03 15:20:40,53] [info] Shutting down WorkflowManagerActor - Timeout = 3600 seconds
[2020-12-03 15:20:40,53] [info] WorkflowManagerActor All workflows finished
[2020-12-03 15:20:40,53] [info] WorkflowManagerActor stopped
[2020-12-03 15:20:40,53] [info] Connection pools shut down
[2020-12-03 15:20:40,53] [info] Shutting down SubWorkflowStoreActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down JobStoreActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down CallCacheWriteActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] SubWorkflowStoreActor stopped
[2020-12-03 15:20:40,53] [info] Shutting down ServiceRegistryActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] Shutting down DockerHashActor - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] JobStoreActor stopped
[2020-12-03 15:20:40,53] [info] Shutting down IoProxy - Timeout = 1800 seconds
[2020-12-03 15:20:40,53] [info] CallCacheWriteActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] WriteMetadataActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] KvWriteActor Shutting down: 0 queued messages to process
[2020-12-03 15:20:40,53] [info] CallCacheWriteActor stopped
[2020-12-03 15:20:40,53] [info] DockerHashActor stopped
[2020-12-03 15:20:40,53] [info] IoProxy stopped
[2020-12-03 15:20:40,53] [info] ServiceRegistryActor stopped
[2020-12-03 15:20:40,54] [info] Database closed
[2020-12-03 15:20:40,54] [info] Stream materializer shut down
Workflow 837d7e7d-ec48-4cfc-8f9e-4736c8bc5a8a transitioned to state Failed
[2020-12-03 15:20:40,57] [info] Automatic shutdown of the async connection
[2020-12-03 15:20:40,57] [info] Gracefully shutdown sentry threads.
[2020-12-03 15:20:40,57] [info] Shutdown finished.

```

# Answer
I updated to the newest version of Cromwell and the error went away!
> 1 votes

---
Tags: gatk, cnv
---


title: Block wise protein imputation

text:
# Question
Title: Block wise protein imputation
I am currently working on a dataset that contains 50 samples (10 samples * 5 blocks). The features of the date set are:

1. The data is perfectly balanced between blocks, with equal treatment representation in each block. Each block contains 2 control samples (CTL) that are pools of all other treatment samples. These are technical replicates that are identical across multiplexes.
2. There is a lot of missing data (~40%),
3. most of the missing data is missing by block, meaning that the protein is present in some blocks, but not others.
4. When looking at the expression of these missing proteins in the other blocks, they tend to be around the median value for the experiment, which suggests to me that these values are missing completely at random.

I could proceed using complete cases, but then I'm missing out on a large portion of data. I am considering using an imputation method, but I'm concerned that (due to the nature of the missing data) that I will just be introducing a block effect into the data.

Any input on potential methods for imputation would be greatly appreciated.

Experimental design:

```
SampleNumber    Multiplex   Treatment
1   A   CTL
2   A   CTL
3   A   TMT1
4   A   TMT1
5   A   TMT2
6   A   TMT2
7   A   TMT3
8   A   TMT3
9   A   TMT4
10  A   TMT4
11  B   CTL
12  B   CTL
13  B   TMT1
14  B   TMT1
15  B   TMT2
16  B   TMT2
17  B   TMT3
18  B   TMT3
19  B   TMT4
20  B   TMT4
21  C   CTL
22  C   CTL
23  C   TMT1
24  C   TMT1
25  C   TMT2
26  C   TMT2
27  C   TMT3
28  C   TMT3
29  C   TMT4
30  C   TMT4
31  D   CTL
32  D   CTL
33  D   TMT1
34  D   TMT1
35  D   TMT2
36  D   TMT2
37  D   TMT3
38  D   TMT3
39  D   TMT4
40  D   TMT4
41  E   CTL
42  E   CTL
43  E   TMT1
44  E   TMT1
45  E   TMT2
46  E   TMT2
47  E   TMT3
48  E   TMT3
49  E   TMT4
50  E   TMT4

```

The data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression.

The controls are technical replicates, but they are pools of all other 40 samples. This is for TMT proteomics, where the pools are used by the MS to identify peptides across the multiplexes, not necessarily for multiplex correction. I've included multiplex (as well as treatment) in the linear model to identify changes associated with the treatment, while removing the effect of the multiplex. I have confirmed (with PCA and PVCA) that the multiplex effect is negligible after including it in the lm. My problem at the moment is the missing data.

I calculated the missing percentage from the data matrix by taking the number of missing values divided by the total number of potential values. he data is in a matrix format, columns = samples and rows = individual proteins, similar to gene expression. To calculate the missing values, I divided the total number of NA's by the total number of values in the matrix (number of samples * number of proteins). Missing data is common in proteomics, but normally this can be attributed to something, most often being close to the limit of detection. But that is not the case here, if a value is missing, it is missing in the whole multiplex, including the pool, where we know it should be found.

The missing values in this case are random, and this is definitely due to the MS. It will select the first ~65k peptides that it identifies to trigger an additional MS event, so if it doesn't see that first peptide, it will be missing from the entire multiplex. There are methods for imputation of data that is missing at random, but I'm concerned that it will just be imputing by multiplex.

# Answer
40% missing data is huge. Missing data analysis is complicated on the underlying distribution. If the data set is periodic then missing data periodicity is needed. Non-periodic data can be solved using regression, I think `SciPy` has an automated method for this, but the volume of missing data makes this approach complex. If you can bring the missing data down to 20% by compartmenatlising the data set that would be the starting point.
> 1 votes

# Answer
This looks like a TMT proteomics experiment. What tools/pipeline did you use to arrive at your quantitative values?

There are a variety of tools that are useful for analyzing this type of data and supply appropriate imputation methods. If you are comfortable programming in R, one I would suggest is the MSstatsTMT R package: http://msstats.org/msstatstmt/

That being said, any ways you can improve the upstream analysis to decrease the amount of missing values will be extremely beneficial.
> 1 votes

---
Tags: statistics, imputation, proteomics
---


# Question
Title: Is there a straightforward way to get mismatches/indels in a BAM file using pysam?
Ideally, I'd like to do something like the following:

```
f = pysam.AlignmentFile(filename, index_filename=index_filename)
read = list(f.fetch('Chromosome', 1, 10000))[0]
read.get_variations()

```

And have that output something like:

```
[
  (1, 'A', 'C'), # An A -> C mismatch at position 1
  (15, 'I', 'CGGG'), # An insertion of the string 'CGGG' starting at position 15 of the read
  (32, 'D', 4), # A deletion of 4 base pairs starting at position 32
]

```

I know that it's possible to figure out all of those values by parsing the CIGAR string and MD tags but I'm wondering if there's a more intuitive interface.

# Answer
> 1 votes
You're looking for `pileup`, which is the htslib (and thus samtools/bcftools) method for finding variants. See this section of the pysam documentation. You may want to peruse about the output as well. The gist is:

```
import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb" )
for pileupcolumn in samfile.pileup("chr1", 100, 120):
    for pileupread in pileupcolumn.pileups:
        if not pileupread.is_del and not pileupread.is_refskip:
            # query position is None if is_del or is_refskip is set.
            ...handle mismatches...

```

---
Tags: sam
---


title: How to obtain desired output?

text:
# Question
Title: How to obtain desired output?
I am working on a project using the following command within `nano`:

```
from Bio import SeqIO
import sys
import re

    fasta_file = (sys.argv[1])
    for myfile in SeqIO.parse(fasta_file, "fasta"):
      if len(myfile) > 250:
       gene_id = myfile.id
       mylist = re.match(r"H149xcV_\w+_\w+_\w+", gene_id)
       print (">"+mylist.group(0))   # edited from list to mylist

```

and its providing with the following outout:

```
>H149xcV_Fge342_r3_h2_d1
>H149xcV_bTr423_r3_h2_d1
>H149xcV_kN893_r3_h2_d1
>H149xcV_DNp021_r3_h2_d1
>H149xcV_JEP3324_r3_h2_d1
>H149xcV_SRt424234_r3_h2_d1

```

How can I change my command so that it provides me with the desired format and only the **UNIQUE** gene id's:

```
>H149xcV_Fge342_r3_h2
>H149xcV_bTr423_r3_h2
>H149xcV_kN893_r3_h2
>H149xcV_DNp021_r3_h2
>H149xcV_JEP3324_r3_h2
>H149xcV_SRt424234_r3_h2

```

# Answer
> 2 votes
its easy replace `match` with `sub` but please stop using `list` as a variable ... mylist is fine.

This might work

```
   mylist = re.sub(r'H149xcV_\w+_\w+_\w+', gene_id)

```

otherwise

```
myregex = re.compile('_\w+\s+.*')
fastaid = myregex.sub('', myfile)

```

OR from @MaximilianPress

```
myregex2 = re.compile('_\w+\n') # or myregex2 = re.compile('_\w+$')
fastaid2 = myregex.sub('\n', myfile) # or fasaid2 = myregex.sub('', myfile)

```

The above will work ... As with all my code I never proof it ..

# Answer
> 0 votes
```
from Bio import SeqIO
import sys
import re

unique = []
fasta_file = (sys.argv[1])
for myfile in SeqIO.parse(fasta_file, "fasta"):
    if len(myfile) > 250:
        gene_id = myfile.id
        if gene.id not in unique:
            unique.append(gene.id)
            mylist = re.match(r"H149xcV_\w+_\w+_\w+", gene_id)
            print (">"+mylist.group(0))

```

let me know if it works I am learning too

---
Tags: python, phylogenetics, fasta, phylogeny
---


title: What is the origin of HIV1?

text:
# Question
Title: What is the origin of HIV1?
What is the origin of HIV1?

This is the alignment of a small part of the POL protein of the HIV1 virus GenBank KU749412.1, with the POL protein of the Visna virus GenBank L06906.1

```
IKKKDSTKWRKLVDFRELNKRTQDFWEVQLGIPHPAGLKKKKSVTVLDVGDAYFSVPLDESFRKYTAFTI--PSTNNETPGIRYQYNVLPQGWKGSPAIFQSSMTKILEPFRAQNPEIVIYQYMDDLYVGSDLEIGQHRAKIEELREHLLRWGFTTPDKKHQKEPPF--------LWMGYELHPDKWTVQPIQLPE--KDSWTVNDIQKLV
|||| | ||| | ||||||| | |  | ||| ||| ||  || || || |||||  || |  | || |    |  ||  |  ||   ||||||| |||  |  | |||       | |    |||| | ||||  |      || |                                | |||| ||  |   |||      | |  ||||
IKKK-SGKWRMLIDFRELNKQTEDLAEAQLGLPHPGGLQRKKHVTILDIGDAYFTIPLYEPYRQYTCFTMLSP--NNLGPCVRYYWKVLPQGWKLSPAVYQFTMQKILRGWIEEHPMIQFGIYMDDIYIGSDL--GL-----EEHRGIVNELASYIAQYGFMLPEDKRQEGYPAKWL-GFELHPEKWKFQKHTLPEITEGPITLNKLQKLV

```

Are these matches,

* By chance or convergence,
* HIV1 and Visna virus share a recent common ancestor or,
* the common ancestor was closer Visna virus from where HIV descended?

Could it be that HIV1 really comes from some kind of experiment involving Visna virus, as I read on the Internet?

The disease caused by the Visna virus is (100% deadly see https://pubmed.ncbi.nlm.nih.gov/16425963/ ). Could there be parallels to sars2?

But look at it: the Visna virus is transmitted by air the first symptoms are like pneumonia (or a brain infection) and otherwise are chronic progressing slowly like SARS-CoV-2 infection resulting in long Covid.

Symptoms are also felt after many days that the virus has come into action, this seems to me not to happen for coronaviruses or for flu viruses but for HIV and Visna virus.

Related links:

1. https://pubmed.ncbi.nlm.nih.gov/8289690/
2. https://pubmed.ncbi.nlm.nih.gov/16425963/
3. (preprint sars2 as retrovirus) https://www.biorxiv.org/content/10.1101/2020.12.12.422516v1
4. https://www.gao.gov/assets/100/91319.html

# Answer
> 6 votes
They are both lentiviruses and share a distant common ancestor.

HIV-1 and HIV-2 are descendents from simian immunodeficiency virus (SIV). The following tree shows the relationships very clearly, from Wertheim and Worobey (2009) Dating the Age of the SIV Lineages That Gave Rise to HIV-1 and HIV-2 PLoS Comp Biol. here. The evolution of HIV-1 is heavily mixed with chimanpzee (cpz) SIV, and there are regions of highly anomalous homology juxtaposed with regions of notably lower homology.

---
Tags: phylogenetics, proteins, phylogeny, covid-19, virus
---


# Question
Title: annotation using ChIPseeker package error
I have differential binding sites object obtain from diffBind (dba.report). I am using the ChIP Seeker package to annotate the peaks but keep getting the following error: Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’

I used diffloop package to add chr to the chromosome column but it did not work. Result in the same error.

Please let me know if you have any idea how to fix it. Highly appreciated Here is my peakset

BD

  seqnames              ranges strand |      Conc Conc_condition Conc_Control      Fold      p-value         FDR
               <Rle>           <IRanges>  <Rle> | <numeric>     <numeric>    <numeric> <numeric>    <numeric>   <numeric>
    [1] NC_000003.12   75668856-75670016      * |  12.27197    12.8243533      0.76175  12.06260 1.43013e-102 1.41852e-97
    [2] NC_000024.10   11320736-11323578      * |  10.58988    11.1419124      0.76175  10.38016 1.15927e-100 5.74930e-96
    [3] NC_000004.12         37431-38634      * |   6.42233     0.0596816      8.06241  -8.00273  6.67673e-98 2.20750e-93
    [4] NC_000017.11   21678164-21686752      * |   9.26507     9.8163416      0.76175   9.05459  3.54723e-95 8.79607e-91
    [5] NC_000022.11   12692637-12694102      * |   9.64945    10.2010195      0.76175   9.43927  2.35838e-87 4.67845e-83
    ...          ...                 ...    ... .       ...           ...          ...       ...          ...         ...
  [226] NC_000004.12 131722954-131744030      * |   11.8156      12.36301     5.337173   7.02584  2.54289e-17 4.78151e-16
  [227] NC_000015.10   68134159-68136364      * |    8.2514       8.80099     0.969455   7.83153  2.28289e-16 3.73532e-15
  [228] NC_000005.10   34189082-34194027      * |   10.4671      11.01445     3.990668   7.02378  1.46745e-15 2.12829e-14
  [229] NC_000018.10       108003-108613      * |   11.2783       5.46414    12.912836  -7.44869  3.68643e-12 3.26794e-11
  [230] NC_000004.12   49709306-49710113      * |   12.0662       6.40350    13.698765  -7.29526  1.23922e-08 6.37694e-0

peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
> peakAnno <- annotatePeak(BD, tssRegion=c(-3000, 3000), TxDb=txdb, annoDb="org.Hs.eg.db")
>> preparing features information...         2020-11-12 01:08:37 PM
>> identifying nearest features...       2020-11-12 01:08:37 PM
>> calculating distance from peak to TSS...  2020-11-12 01:08:38 PM
>> assigning genomic annotation...       2020-11-12 01:08:38 PM
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘NSBS’ for signature ‘"SortedByQueryHits"’```
</code></pre>

# Answer
that was solved by changing seq level:

<pre><code>newnames <- paste0(c("1","2","3","4", "5","6","7", "8","9","10","11","12","13","14","15","16","17","18","19","20","21","22","X","Y","MT"))
names(newnames) <- paste0(c("NC_000001.11","NC_000002.12","NC_000003.12","NC_000004.12", "NC_000005.10","NC_000006.12","NC_000007.14", "NC_000008.11","NC_000009.12","NC_000010.11","NC_000011.10","NC_000012.12","NC_000013.11","NC_000014.9","NC_000015.10","NC_000016.10","NC_000017.11","NC_000018.10","NC_000019.10","NC_000020.11","NC_000021.9","NC_000022.11","NC_000023.11","NC_000024.10","NC_012920.1"))

peaks_with_new_seq_level <- renameSeqlevels(peaks_with_old_seq_level,newnames)```

> 2 votes

---
Tags: differential-expression, chip-seq, sequence-annotation, peak-calling
---


# Question
Title: the model matrix is not full rank : This is a classic question which a biologist face without clear understanding of the model design
Saw this answer biostar. Tried to make my metadata as such but still the

```
"Error in checkFullRank(modelMatrix) : "

```

This is my coldata

```
dput(coldata)
structure(list(Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L,
6L, 6L, 6L, 6L), .Label = c("HSC", "Blast", "CMP", "GMP", "LSC",
"Mono"), class = "factor"), Mutation = structure(c(10L, 4L, 3L,
5L, 2L, 3L, 9L, 11L, 1L, 1L, 7L, 7L, 10L, 7L, 8L, 2L, 10L, 10L,
10L, 10L, 2L, 2L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 5L, 2L,
9L, 1L, 1L, 8L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("DNMT3A-R882H",
"FLT3-ITD", "IDH2-R140Q", "KRAS-G13D", "MAU2-Q98H", "No", "NPM1-L287ins(TCTG)",
"NRAS-G13D", "STAG2-R614", "TET2-E1357", "TET2-R550"), class = "factor"),
    Satus = structure(c(2L, 5L, 4L, 6L, 3L, 4L, 9L, 10L, 1L,
    1L, 7L, 7L, 2L, 7L, 8L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 7L, 16L,
    16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
    16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 12L, 12L, 13L,
    12L, 9L, 11L, 11L, 14L, 15L, 15L, 15L, 15L, 15L, 15L), .Label = c("Blast-DNMT3A-R882H",
    "Blast-E1357", "Blast-FLT-ITD", "Blast-IDH2-R140Q", "Blast-KRAS-G13D",
    "Blast-MAU2-Q98H", "Blast-NPM1-L287ins(TCTG)", "Blast-NRAS-G13D",
    "Blast-STAG2-R614", "Blast-TET2-R550", "LSC-FLT-DNMT3A-R882H",
    "LSC-FLT-ITD", "LSC-MAU2-Q98H", "LSC-NRAS-G13D", "Mature-Normal",
    "Progenitor-Normal", "Stem-Normal"), class = "factor")), row.names = c("Blast1",
"Blast10", "Blast11", "Blast12", "Blast13", "Blast14", "Blast15",
"Blast16", "Blast17", "Blast18", "Blast19", "Blast20", "Blast2",
"Blast21", "Blast22", "Blast23", "Blast3", "Blast4", "Blast5",
"Blast6", "Blast7", "Blast8", "Blast9", "CMP1", "CMP2", "CMP3",
"CMP4", "CMP5", "CMP6", "CMP7", "CMP8", "GMP1", "GMP2", "GMP3",
"GMP4", "GMP5", "GMP6", "GMP7", "HSC1", "HSC2", "HSC3", "HSC4",
"HSC5", "HSC6", "HSC7", "LSC1", "LSC2", "LSC3", "LSC4", "LSC5",
"LSC6", "LSC7", "LSC8", "Mono1", "Mono2", "Mono3", "Mono4", "Mono5",
"Mono6"), class = "data.frame")

dds <- DESeqDataSetFromMatrix(countData=df2, colData=coldata, design= ~ Group + Mutation)

```

This is the design I'm trying to incorporate. What i get is my design is nested in one of the column by looking at numerous example.

What i would like to do is I'm trying to run a LRT test on this dataset which works fine if i give "`Group`" as my design and then i run this

```
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

```

which gives me diffrenitaly expressed across multiple groups. Now there are two question

1. Is it possible to include interaction such as in my case `Group+Mutation` and do the lrt test?
2. Or the above is conceptually wrong to do ?

# Answer
```
> table(doldat$Group, coldat$Mutation)

        DNMT3A-R882H FLT3-ITD IDH2-R140Q KRAS-G13D MAU2-Q98H No NPM1-L287ins(TCTG) NRAS-G13D STAG2-R614 TET2-E1357 TET2-R550
  HSC              0        0          0         0         0  7                  0         0          0          0         0
  Blast            2        4          2         1         1  0                  4         1          1          6         1
  CMP              0        0          0         0         0  8                  0         0          0          0         0
  GMP              0        0          0         0         0  7                  0         0          0          0         0
  LSC              2        3          0         0         1  0                  0         1          1          0         0
  Mono             0        0          0         0         0  6                  0         0          0          0         0

```

Look at this table, you see by eye how some mutations are nested with Group. If you want to include it then you will need to subset the experiment into multiple separate runs, depending on the question you want to answer.
> 3 votes

---
Tags: r, deseq2, design
---


# Question
Title: Biology behind PCA analysis based on SNP
What is the biology behind SNP based PCA analysis to study population structure? I am reading some articles where they used PCA analysis to compare isolates of drosophila that is collected from a different location. AND by PCA they found that SNP is different between all populations. What does that mean? (what it tell you about population structure/genetics?)

# Answer
First of all, PCA is a technique for dimension reduction. Basically, the goal is to compare tens of thousands of SNPs in Drosophila. Now if you only have 2 SNPs, you can plot them on a 2D scatter plot. If you have 3 SNPs, you may try a 3D plot. But now imagine you have 30,000 SNPs, but you CANNOT plot a 30000-dimensional plot. To visualize this high dimensional data, what we can do is to perform dimensional reduction like PCA. PCA tries to find a set of orthogonal coordinations that explains most of the variation in the data (if there no variation, there is no information contained in the data, which essentially means there is no data). The idea is that PC1 carries most variation can be explained, and PC2 carries the second most... For lower PCs like PC50 or PC60, they probably only carry noise in the data. Therefore, the higher PCs (PC1, PC2 and so on) effectively summarizes the useful information in the data. So you can visualize the "structure" of the data in a 2D PCA plot.

By looking at the distance between points on a PCA plot, you can tell how similar the two data points are. But if you see two populations that are perfectly separated on PCA plot, it does not mean that the 2 population differ completely at every SNP, because PCA is a summarization of all SNP included.
> 1 votes

---
Tags: snp, pca
---


# Question
Title: What unit is the length of a protein measured in?
Sorry if this question is obvious - I'm not familiar with proteomics so I'm not sure of the standard units.

In this paper, it talks about using the length of a protein to calculate the NSAF value. However, it doesn't specify the units. Would it be peptides, amino acids or nucleic acids? Or something else?

# Answer
As mentioned by others, it's the number of amino acids. Alternative measurements could be in Angstroms or similar, but that would require at least an approximate structure.
> 1 votes

---
Tags: proteins, proteomics
---


# Question
Title: How to select specific lines and specific columns
I have a data frame that contain lot of columns and lines, like this for example

```
Column1   Column2
5            T1
1            T2
1            T3
1            T1
2            T1
1            T1
2            T1
1            T1

```

I want to extract just lines that have the value 1 in column 1 and Value T1 in Column 2

So the result I want is

```
Column1   Column2
1            T1
1            T1
1            T1

```

Can you tell me please how I can do that

# Answer
You can simply use subset function in R and use the code given below,

```
    new_dataframe <-
    subset(name_of_your_dataframe,Column1 == 1 &&
    Column2 ==  "T1")

```

Please keep in mind that R is case sensitive and you should be careful while using the variable names.
> 2 votes

---
Tags: r, statistics, database
---


# Question
Title: ggplot2 geom_bar() , do not merge y-values shared same x-ids
I am plotting a barplot using `ggplot` \+ `geom_bar()` function in R.

I have some the recoreds which share the same x-axis identifiers.

`geom_bar()` did merge all values that shared the same x-axis identifiers.

Like the following example.

```
library("ggplot2")

df <- data.frame(x = c("a","a","b","c","d","e"),
                 y = c(5,5,5,5,5,5),
                 stringsAsFactors = FALSE)

ggplot(df, aes(x, y)) + geom_bar(stat = "identity")

```

The erroneous graph is displayed below.

The highest bar means the sum of multiple y-values shared same x-axis names.

Anyone else could help me to solve this problem? Thanks.

# Answer
> 1 votes
As you seem to have surmised, the graph is not erroneous, but correct. You have a total of 10 for `a` and that's what it's showing. If you want it to ignore duplicate X-axis values then you need to tell ggplot which datapoints it should actually be using. `barplot()` is assuming that your duplicate `a` values are intentional, ggplot assumes they should be grouped. These are both reasonable assumptions and it's your responsibility to determine beforehand which is correct.

---
Tags: r, ggplot2
---


title: How can I convert Protein ID to PDB ID?

text:
# Question
Title: How can I convert Protein ID to PDB ID?
I have a list of Protein ID such as Q9S9U6, A2XLL2, Q06402. What name is this ID? UniProtKB ID? And how can I get PDB ID from this list? Is its PDB ID only one?

I'm using Python.

Get a list of PDB ids from a list of sequence protein number?

That method of this link gives me some IDs but not only one and usually empty.

# Answer
Get their sequences and then search the pdb (by his API) for the obtained sequences

https://search.rcsb.org/#search-services

...

sequence This service employs the MMseqs2 software and performs fast sequence matching searches (BLAST-like) based on a user-provided FASTA sequence (with E-value or % Identity cutoffs). Following search targets are available: pdb\_protein\_sequence: all current protein sequences in PDB...
> 1 votes

# Answer
I think this should help you : https://biit.cs.ut.ee/gprofiler/convert
> 0 votes

---
Tags: python, proteins
---


# Question
Title: Which one is a more convenient assembly?
I have developed a software for de novo genome assembly. Its performance varies gradually according to how much data you employ. At initial stages it often produces contigs that look like that when aligned to a reference genome:

**Initial-stage s-aligner assembly of SRR8357422** Genome fraction:98,8% - Largest Alignment: 61.544 - NGA50: 59.409 - LGA50: 2

If I add more data, the contigs end up being larger, but that obviously takes more time. My question is, how does this result compare to, for example, one like that?

**SPAdes assembly of SRR8357422** Genome fraction:93,46% - Largest Alignment: 128.391 - NGA50: 128.391 - LGA50: 1

Note how you can construct a contig covering almost the full reference genome assembling green contigs in the first case. Therefore, a blast search of genes will likely let you find more genes given that it covers more genome fraction.

Apart from determining what option is more convenient, I would also like to know how much more convenient is to have both results. Or in which cases it's more convenient.

Last (but not least!), is there a way to generate automatically larger contigs from these overlapping contigs in assembly 1? I have tried using miniasm with no success. Even minimap2 seems not to find overlaps well enough. I am developing my own algorithm also for that but would like to know if there is another way.

\[EDITED on 17/12 to add another example with real data for both assemblies\]

# Answer
I asked that on Reddit. I didn't get an answer to all my questions but I got an answer for how to assemble the overlapping contigs.

If I use Flye I can assemble these and get larger contigs, better NG50, and higher genome fraction.

**s-aligner + Flye assembly of SRR8357422** Genome fraction:95,95% - Largest Alignment: 134.103 - NGA50: 134.103 - LGA50: 1

Given that the data in the first assembly apart from being more useful itself in my opinion, is objectively more suitable to generate a better result in all senses with a simple additional step that takes a few seconds, my only possible conclusion is that it is overall a better result and therefore this sort of result will often be preferable.

Of course, there are other aspects to consider, like error rate, that are not considered here as it wasn't part of the question (mainly because I think Quast is not adequate for calculating it). But BTW, mismatches per 100kbp are also lower.

\[EDIT\] I want to add that this is not granted that in all cases Flye will generate a clearly better assembly.
> 0 votes

---
Tags: assembly, reference-genome, genome-sequencing, scaffold
---


# Question
Title: How to merge back the processes by sample after splitting by chromosome in Nextflow?
I have a pipeline that has an input list of samples. I would like to split those samples by chromosomes, process each chromosome apart, and then merge them back by sample name.

Here's an example:

```
samplesList = Channel.from(["sample1", "sample2"])

// A process that splits the samples by chromosome as an example
// In my actual case, I have a VCF that is split by `snpSift Split` command
// and generates files like [sample]_[chr].vcf
// In the example below, I only gave sample1
// But take into consideration that I have multiple samples

process bar {
   echo true

  output:
    file("*.tsv") into bar_output

  script:
  """
  echo 1 > sample1_chr1.tsv
  echo 2 > sample1_chr2.tsv
  echo 3 > sample1_chr3.tsv
  echo 4 > sample1_chr4.tsv
  echo 5 > sample1_chr5.tsv
  """
}

// Process each chromosome apart for each sample
// This means the output of the process foo
// is only ONE chromosome per sample
process foo {
  input:
    each x from bar_output  // to treat each file separately from the channel
  output:
    file("*.tsv") into xasxa
  script:
    """
    echo \\\$(cat $x)\\\$(cat $x) > sample1_chr\\\$(cat $x)_2.tsv
    """
}

// How should I write the next process to group back them together by sample name?
// I do not wish to do a collect() on channel foo because it will force
// all the processes of foo to finish first.
// In this example, foo would emit 5 different channels of chromosome for sample1
// This process should merge sample1 only when all chromosomes are processed
// independently of other samples.
process mergeChromosomesBySamples {
   input:
     ???
}

```

As the example above, the process `mergeChromosomesBySamples` should return something similar/close to

```
["sample1", ["/path/to/workdir/sample1_chr1.tsv",
             "/path/to/workdir/sample1_chr2.tsv",
             "/path/to/workdir/sample1_chr3.tsv",
             "/path/to/workdir/sample1_chr4.tsv",
             "/path/to/workdir/sample1_chr5.tsv"]
]

```

---

EDIT1: more explication for the process `mergeChromosomesBySamples` is given in the comments and as follows ...

The channel `xasxa` processes each chromosome apart, this means that the output of the channel `xasxa` contains only one chromosome of the sample like `"/path/to/workdir/sample1_chr1.tsv"` only. My goal is then to merge the processes out `xasxa` for each sample together.

I do not wish to use `collect()` in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample.

Where each sample has its own VCF generated and I process them independently.

EDIT2: Here's a graph of what I am trying to achieve:

If I use `.collect()`, then all the processing on the chromosomes of all the samples has to finish first before being able to merge them, which is something I do not want.

I appreciate your help.

# Answer
> 2 votes
I was able to find a solution that answers my problem finally.

So, here are two different methods below, one which does exactly what I want while the other has to wait for all the previous processes to finish:

**First method**: waiting for previous process to finish using `.groupTuple()` like what @Pallie mentionned. This method is simpler than the second method. The problem is that it has to wait for the previous processes to finish. I want my samples to be treated in parallel. If I missed something with the logic of the method, please feel free to edit.

```
#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
         echo 2 > ${x}_chr2.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr*vcf") into bar_out
      script:
      """
         echo ${x}
         filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x) from bar_out.groupTuple()

      script:
      """
         echo $sample
         echo ${x}
      """
}

```

will output:  As you can see above, the process `mergefoobar` is still waiting for all the processes of `bar` to finish.

**Second method**: does not wait for the previous processes to finish, uses `.join()` for each chromosome... Yes, the downside is that you would have to write the `join()` for all the chromosomes. Additionally, if you use `snpSift Split`, one sample may generate 3 files out of 5 (because only 3 chromosomes had variants listed in them), you would need to create the other empty VCF files because the output excepts them.

```
#!/usr/bin/env nextflow

// my samples
samplesList = Channel.from(["sample1", "sample2"])

// the process that is supposed to split the sample by chromosomes
process foo {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         val(x) from samplesList
      output:
         file("*.vcf") into foo_out
         val(x) into output_sample
      script:
      """
         echo 1 > ${x}_chr1.vcf
     echo 2 > ${x}_chr2.vcf
         echo X > ${x}_chrX.vcf
     echo Y > ${x}_chrY.vcf
      """
}

//Combine the samples by chromosome and filter out only the samples that match the files output
foo_out.flatten().combine(output_sample).filter{ it[0].toString().contains(it[1].toString()) }.set{maybework}

// the process that is supposed to treat each chromosome apart
process bar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
         set file(y), val(x) from maybework
      output:
         set val(x), file("*chr1*vcf") optional true into bar_out1
         set val(x), file("*chr2*vcf") optional true into bar_out2
         set val(x), file("*chrX*vcf") optional true into bar_outX
         set val(x), file("*chrY*vcf") optional true into bar_outY
      script:
      """
         echo ${x}
     filename=\\\$(basename -- ${y})
         chromosome=\\\$(echo "\\\${filename%.*}" | cut -f2 -d'_')
         cat ${y} > ${x}_"\\\$chromosome"_"\\\$chromosome".vcf
         if [ ${x} == "sample2" ]; then sleep 20; fi
      """
}

// the process that merges back the chromosomes by sample
process mergefoobar {
      echo true
      cpus 1
      memory '1 MB'
      executor 'local'

      input:
      set val(sample), file(x), file(y), file(i), file(t) from bar_out1.join(bar_out2).join(bar_outX).join(bar_outY)

      script:
      """
         echo $sample
     echo ${x}
         echo ${y}
     echo ${i}
         echo ${t}
      """
}

```

Above, you see that the process mergefoobar has launched even if the previous bar has not finished yet.

I hope this also helps anyone in the future who might have this issue.

# Answer
> 5 votes
Use the groupBy operator. If the content of channel `xasxa` is:

`["/path/to/workdir/sample1_chr1.tsv", "/path/to/workdir/sample1_chr2.tsv", "/path/to/workdir/sample2_chr1.tsv", "/path/to/workdir/sample2_chr2.tsv`

You could probably do:

`xasxa.groupBy{ it -> it.fileName.toString().split("_")[0] }`

The result would be a channel with:

`[sample1:[/path/to/workdir/sample1_chr1.tsv, /path/to/workdir/sample1_chr2.tsv, sample2:[/path/to/workdir/sample2_chr1.tsv, /path/to/workdir/sample2_chr1.tsv]]`

This breaks if you have an underscore in your samplename before `_chr*`. More robust would be to initialise all channels with a sample value, so that each channel is in the form of \[samplename, chromosome\]. Then all you need is a `groupTuple` on your channel.

---
Tags: nextflow
---


# Question
Title: Functional Annotation vs Functional Enrichment--which one should I use for Network Analysis?
I am currently working on a network pharmacology project where I am supposed to find out the possible diseases our drug might have an effect on. I am not a bioinformatics student, so my knowledge in the subject matter is very limited. I have all the targets and I am supposed to get the pathways associated with the targets. I used `STRING` for the purpose.

* I got both `KEGG` functional annotation and enrichment for my list of targets.
* I have to build a two-mode network

**Question** I am not very clear about

* the difference between enrichment and functional annotation and
* which data should be used for the next step.

# Answer
> 1 votes
The answer to your main doubt about "annotation vs. enrichment" is that you need both, but you already have one annotation.

You can only perform a functional enrichment analysis after you have a functional annotation. In your question, you state that you have KEGG annotation for your targets. That is it, knowing which target "belongs" to each KEGG module is an example of functional annotation of your targets. Several tools can annotate genes, proteins, etc. Knowing which protein is related to each drug or disease is also a functional annotation.

Now, knowing which targets are significantly important can be achieved by functional enrichment analysis.

For the network, both information can be used together. For example, you could link all elements to their functional annotation and use the available enrichment score as a quantitative parameter of importance for calculating correlations.

---
Tags: proteins, networks, go-enrichment, functional-annotation, kegg
---


# Question
Title: How to count RNA sequencereads using custom made python scripts?
I am trying to do RNA seq analysis and my goal is to filter gene counts less than 5 using custom made python scripts. The code chunk goes as follows

```
#!/usr/bin/python
#Load relevant modules
import sys

#Read in the HTseq gene count output files for each of the 8 samples
inputlist = [sys.argv[1], sys.argv[2], sys.argv[3], sys.argv[4], sys.argv[5], sys.argv[6], sys.argv[7], sys.argv[8]]
##Determine the genes with sufficient coverage
#Set up a variable to hold genes with sufficient coverage in any of the 8 samples
genelist = []

#For each HTseq file
for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')
    #For each line
    for line in inputfile:
        #Determine the number of reads/transcripts for this gene
        splitline = line.strip().split('\t')
        #Does it pass the threshold?
        if int(splitline[1]) >= 5:
            #If so, add it to list
            if not splitline[0] in genelist:
                genelist.append(splitline[0])

    #Close the file
    inputfile.close()

#Write out the list of sufficiently expressed genes
outputlist = open(sys.argv[9], 'w')

for gene in genelist:
    outputlist.write(gene + '\n')

outputlist.close()
##Filter each of the HTseq files for only sufficiently expressed genes across our samples
#For each HTseq file

for i in range(len(inputlist)):
    #Open it
    inputfile = open(inputlist[i], 'r')

#Open an output 'filtered' file
outputfile = open(sys.argv[i+10], 'w')

#For each line
for line in inputfile:
    #Determine the gene
    splitline = line.strip().split('\t')

    #Is it in our list?
    if splitline[0] in genelist:
        #If so, write out the original line
        outputfile.write(line)

#Close the files
inputfile.close()
outputfile.close()

```

`python filtercounts.py 'wth1count.txt' 'wth2count.txt' 'wth3count.txt' 'wth4count.txt' '99h1count.txt' '99h2count.txt' '99h3count.txt' '99h4count.txt' 'genelist.txt' 'wth1filtered.txt' 'wth2filtered.txt' 'wth3filtered.txt' 'wth4filtered.txt' '99h1filtered.txt' '99h2filtered.txt' '99h3filtered.txt' '99h4filtered.txt'`

The program works fine but at the end of the code execution I am not getting the last two output files '99h3filtered.txt' '99h4filtered.txt'.

# Answer
What software are you using for the next step? DESeq in R will do this filtering for you in one line.

I think a more typical way to do this would be to make your entire gene count matrix, and then remove rows where the max value for that row is \< 5. Your way, you have to go through each file twice; once to know which genes to keep, and then again to actually remove lines.
> 1 votes

---
Tags: genomics, rna
---


# Question
Title: How to annotate optimally a fungal genome without RNA-seq evidence?
**Genome information：**

1. ~50M nt
2. 2300+ contigs
3. No pre-trained parameters in Augustus
4. There are several well annotated RefSeq genomes of this genus.
5. There are several RNA-seq data of this genus, but not in this species.

**My current strategy:**

EST evidence:

1. I selected many RNA-seq data that include as many species as possible, and de novo assembled these data using `Trinity`.
2. Only the longest transcripts were saved.
3. All transcript files were catenated into a `Total_transcript.fasta` file, and `cd-hit-est` was used to reduce data size with options `-c 0.8 -n 5 -s 0.8`.

Protein evidence

4. All fungal proteins from SwissProt was download and reduced using `cd-hit` with options `-c 0.8 -n 5 -s 0.8`.
5. Proteinsets from RefSeq genomes of this genus were download, catenated, reduced using same methods.

I planed to annotate genome using `Maker2` with above Evidence.

I have a doubt about my strategy: Maker2 calls Augustus to ab-initio annotation, but no pre-trained model. How should I do to solve this?

Anyone else could share your valuable suggestions? Thanks anyway.

# Answer
The Maker documentation does include some instructions for training *ab initio* gene predictors, but it assumes an abundant EST database is available. (Assumptions about what kinds of sequences will be available for a draft genome assembly have changed drastically in the last decade.)

It may be worth exploring an iterative approach in any case. You can do a first-pass annotation with Maker using the evidence you discussed, along with *ab initio* predictions from Augustus (and I'd recommend SNAP as well). Each gene model will be scored by Maker with the *annotation edit distance (AED)*, representing how well the gene model agrees with the evidence (0.0 is best). From this first pass, if you can identify a few hundred reliable gene models with good AED scores (or maybe even as many as one or two thousand), that should be plenty of training data. You could use these gene models to train a new model for Augustus, using Maker's instructions or Augustus' own instruction. (In my experience, it was quite a bit of work in either case.)
> 1 votes

---
Tags: sequence-annotation
---


# Question
Title: How to write multi-sequence alignment data frame as a fasta format fast
I have a data frame where the row names are the species names, and each column can either be an Amino-Acid character or "-".

I wish to write it in a fasta format. Simple example of the format is:

```
>UniRef90_A0A1V4QIJ3
------------------------------------------------------------
--------------------------MDISTFSP-----PVAG-----------VE----
------ITAF-------------V------------------------------------
LN--------G-------------------------------------------------
----EKIPA-----------------------------------------IVLSKD----
---LNPFL------------------Q-ELLDQEPAC-----------------------
------------EHDIGTCYG------------------------DA-------------
--------------------------------------------------PC-----IS-
-----------------------GPDAWNHSGFQC----Q---------VC---------
------------GLVKSWNL-------------------VGDFVVYRQE-----------
>UniRef90_A0A0L9V5F8
------------------------------------------------------------
RFHTLFRNEYGHLRVLQRFDQRSKQIQNLENYRLVEFKSKPNTLLLPHHADADFLLVVLN
GRALLTLVNP-------------D---------------------------------G--
RDSYILEQGHA-------------------------------------------------
----QKIPAGTIFFLVNPNDNENLRI------------------------IKIATPINNP
HRFQDFFL-------------SSTEAQQSYLQG--FS-----------------------
----------------KNVLEASFDSEFKEINRVLFGEEGQQQQGEESQQEGVIVELERE
QIRELIKHA-----------------------------------------KSSSRRSLS-
-----------------------SQ-----DEPFNLRNRKP--------IY--SN-----
------------KF-GRWYEITPEKNP------------Q--------------------

```

I can do it using the following for loop:

```
for (n in row.names(df)){
    ss = do.call(paste, c(as.list(mu_df[1, ]), sep = ""))
    write.fasta(ss,
                  names=n,
                  as.string=TRUE,
                  open="a",
                  new_file_path)
}

```

Is there a faster way of writing it? I have to write millions of those files and I wish to optimize it.

In order to create the data frame:

```
structure(list(V1 = c("-", "-"), V2 = c("-", "-"), V3 = c("-",
"-"), V4 = c("-", "-"), V5 = c("-", "-"), V6 = c("-", "-"), V7 = c("-",
"-"), V8 = c("-", "-"), V9 = c("-", "-"), V10 = c("-", "-"),
    V11 = c("-", "-"), V12 = c("-", "-"), V13 = c("-", "-"),
    V14 = c("-", "-"), V15 = c("-", "-"), V16 = c("-", "-"),
    V17 = c("-", "-"), V18 = c("-", "-"), V19 = c("-", "-"),
    V20 = c("-", "-"), V21 = c("-", "-"), V22 = c("-", "-"),
    V23 = c("-", "-"), V24 = c("-", "-"), V25 = c("-", "-"),
    V26 = c("-", "-"), V27 = c("-", "-"), V28 = c("-", "-"),
    V29 = c("-", "-"), V30 = c("-", "-"), V31 = c("-", "-"),
    V32 = c("-", "-"), V33 = c("-", "-"), V34 = c("-", "-"),
    V35 = c("-", "-"), V36 = c("-", "-"), V37 = c("-", "-"),
    V38 = c("-", "-"), V39 = c("-", "-"), V40 = c("-", "-"),
    V41 = c("-", "-"), V42 = c("-", "-"), V43 = c("-", "-"),
    V44 = c("-", "-"), V45 = c("-", "-"), V46 = c("-", "-"),
    V47 = c("-", "-"), V48 = c("-", "-"), V49 = c("-", "-"),
    V50 = c("-", "-"), V51 = c("-", "-"), V52 = c("-", "-"),
    V53 = c("-", "-"), V54 = c("-", "-"), V55 = c("-", "-"),
    V56 = c("-", "-"), V57 = c("-", "-"), V58 = c("-", "-"),
    V59 = c("-", "-"), V60 = c("-", "-"), V61 = c("-", "R"),
    V62 = c("-", "F"), V63 = c("-", "H"), V64 = c("-", "T"),
    V65 = c("-", "L"), V66 = c("-", "F"), V67 = c("-", "R"),
    V68 = c("-", "N"), V69 = c("-", "E"), V70 = c("-", "Y"),
    V71 = c("-", "G"), V72 = c("-", "H"), V73 = c("-", "L"),
    V74 = c("-", "R"), V75 = c("-", "V"), V76 = c("-", "L"),
    V77 = c("-", "Q"), V78 = c("-", "R"), V79 = c("-", "F"),
    V80 = c("-", "D"), V81 = c("-", "Q"), V82 = c("-", "R"),
    V83 = c("-", "S"), V84 = c("-", "K"), V85 = c("-", "Q"),
    V86 = c("-", "I"), V87 = c("M", "Q"), V88 = c("D", "N"),
    V89 = c("I", "L"), V90 = c("S", "E"), V91 = c("T", "N"),
    V92 = c("F", "Y"), V93 = c("S", "R"), V94 = c("P", "L"),
    V95 = c("-", "V"), V96 = c("-", "E"), V97 = c("-", "F"),
    V98 = c("-", "K"), V99 = c("-", "S"), V100 = c("P", "K"),
    V101 = c("V", "P"), V102 = c("A", "N"), V103 = c("G", "T"
    ), V104 = c("-", "L"), V105 = c("-", "L"), V106 = c("-",
    "L"), V107 = c("-", "P"), V108 = c("-", "H"), V109 = c("-",
    "H"), V110 = c("-", "A"), V111 = c("-", "D"), V112 = c("-",
    "A"), V113 = c("-", "D"), V114 = c("-", "F"), V115 = c("V",
    "L"), V116 = c("E", "L"), V117 = c("-", "V"), V118 = c("-",
    "V"), V119 = c("-", "L"), V120 = c("-", "N"), V121 = c("-",
    "G"), V122 = c("-", "R"), V123 = c("-", "A"), V124 = c("-",
    "L"), V125 = c("-", "L"), V126 = c("-", "T"), V127 = c("I",
    "L"), V128 = c("T", "V"), V129 = c("A", "N"), V130 = c("F",
    "P"), V131 = c("-", "-"), V132 = c("-", "-"), V133 = c("-",
    "-"), V134 = c("-", "-"), V135 = c("-", "-"), V136 = c("-",
    "-"), V137 = c("-", "-"), V138 = c("-", "-"), V139 = c("-",
    "-"), V140 = c("-", "-"), V141 = c("-", "-"), V142 = c("-",
    "-"), V143 = c("-", "-"), V144 = c("V", "D"), V145 = c("-",
    "-"), V146 = c("-", "-"), V147 = c("-", "-"), V148 = c("-",
    "-"), V149 = c("-", "-"), V150 = c("-", "-"), V151 = c("-",
    "-"), V152 = c("-", "-"), V153 = c("-", "-"), V154 = c("-",
    "-"), V155 = c("-", "-"), V156 = c("-", "-"), V157 = c("-",
    "-"), V158 = c("-", "-"), V159 = c("-", "-"), V160 = c("-",
    "-"), V161 = c("-", "-"), V162 = c("-", "-"), V163 = c("-",
    "-"), V164 = c("-", "-"), V165 = c("-", "-"), V166 = c("-",
    "-"), V167 = c("-", "-"), V168 = c("-", "-"), V169 = c("-",
    "-"), V170 = c("-", "-"), V171 = c("-", "-"), V172 = c("-",
    "-"), V173 = c("-", "-"), V174 = c("-", "-"), V175 = c("-",
    "-"), V176 = c("-", "-"), V177 = c("-", "-"), V178 = c("-",
    "G"), V179 = c("-", "-"), V180 = c("-", "-"), V181 = c("L",
    "R"), V182 = c("N", "D"), V183 = c("-", "S"), V184 = c("-",
    "Y"), V185 = c("-", "I"), V186 = c("-", "L"), V187 = c("-",
    "E"), V188 = c("-", "Q"), V189 = c("-", "G"), V190 = c("-",
    "H"), V191 = c("G", "A"), V192 = c("-", "-"), V193 = c("-",
    "-"), V194 = c("-", "-"), V195 = c("-", "-"), V196 = c("-",
    "-"), V197 = c("-", "-"), V198 = c("-", "-"), V199 = c("-",
    "-"), V200 = c("-", "-"), V201 = c("-", "-"), V202 = c("-",
    "-"), V203 = c("-", "-"), V204 = c("-", "-"), V205 = c("-",
    "-"), V206 = c("-", "-"), V207 = c("-", "-"), V208 = c("-",
    "-"), V209 = c("-", "-"), V210 = c("-", "-"), V211 = c("-",
    "-"), V212 = c("-", "-"), V213 = c("-", "-"), V214 = c("-",
    "-"), V215 = c("-", "-"), V216 = c("-", "-"), V217 = c("-",
    "-"), V218 = c("-", "-"), V219 = c("-", "-"), V220 = c("-",
    "-"), V221 = c("-", "-"), V222 = c("-", "-"), V223 = c("-",
    "-"), V224 = c("-", "-"), V225 = c("-", "-"), V226 = c("-",
    "-"), V227 = c("-", "-"), V228 = c("-", "-"), V229 = c("-",
    "-"), V230 = c("-", "-"), V231 = c("-", "-"), V232 = c("-",
    "-"), V233 = c("-", "-"), V234 = c("-", "-"), V235 = c("-",
    "-"), V236 = c("-", "-"), V237 = c("-", "-"), V238 = c("-",
    "-"), V239 = c("-", "-"), V240 = c("-", "-"), V241 = c("-",
    "-"), V242 = c("-", "-"), V243 = c("-", "-"), V244 = c("-",
    "-"), V245 = c("E", "Q"), V246 = c("K", "K"), V247 = c("I",
    "I"), V248 = c("P", "P"), V249 = c("A", "A"), V250 = c("-",
    "G"), V251 = c("-", "T"), V252 = c("-", "I"), V253 = c("-",
    "F"), V254 = c("-", "F"), V255 = c("-", "L"), V256 = c("-",
    "V"), V257 = c("-", "N"), V258 = c("-", "P"), V259 = c("-",
    "N"), V260 = c("-", "D"), V261 = c("-", "N"), V262 = c("-",
    "E"), V263 = c("-", "N"), V264 = c("-", "L"), V265 = c("-",
    "R"), V266 = c("-", "I"), V267 = c("-", "-"), V268 = c("-",
    "-"), V269 = c("-", "-"), V270 = c("-", "-"), V271 = c("-",
    "-"), V272 = c("-", "-"), V273 = c("-", "-"), V274 = c("-",
    "-"), V275 = c("-", "-"), V276 = c("-", "-"), V277 = c("-",
    "-"), V278 = c("-", "-"), V279 = c("-", "-"), V280 = c("-",
    "-"), V281 = c("-", "-"), V282 = c("-", "-"), V283 = c("-",
    "-"), V284 = c("-", "-"), V285 = c("-", "-"), V286 = c("-",
    "-"), V287 = c("-", "-"), V288 = c("-", "-"), V289 = c("-",
    "-"), V290 = c("-", "-"), V291 = c("I", "I"), V292 = c("V",
    "K"), V293 = c("L", "I"), V294 = c("S", "A"), V295 = c("K",
    "T"), V296 = c("D", "P"), V297 = c("-", "I"), V298 = c("-",
    "N"), V299 = c("-", "N"), V300 = c("-", "P"), V301 = c("-",
    "H"), V302 = c("-", "R"), V303 = c("-", "F"), V304 = c("L",
    "Q"), V305 = c("N", "D"), V306 = c("P", "F"), V307 = c("F",
    "F"), V308 = c("L", "L"), V309 = c("-", "-"), V310 = c("-",
    "-"), V311 = c("-", "-"), V312 = c("-", "-"), V313 = c("-",
    "-"), V314 = c("-", "-"), V315 = c("-", "-"), V316 = c("-",
    "-"), V317 = c("-", "-"), V318 = c("-", "-"), V319 = c("-",
    "-"), V320 = c("-", "-"), V321 = c("-", "-"), V322 = c("-",
    "S"), V323 = c("-", "S"), V324 = c("-", "T"), V325 = c("-",
    "E"), V326 = c("-", "A"), V327 = c("Q", "Q"), V328 = c("-",
    "Q"), V329 = c("E", "S"), V330 = c("L", "Y"), V331 = c("L",
    "L"), V332 = c("D", "Q"), V333 = c("Q", "G"), V334 = c("E",
    "-"), V335 = c("P", "-"), V336 = c("A", "F"), V337 = c("C",
    "S"), V338 = c("-", "-"), V339 = c("-", "-"), V340 = c("-",
    "-"), V341 = c("-", "-"), V342 = c("-", "-"), V343 = c("-",
    "-"), V344 = c("-", "-"), V345 = c("-", "-"), V346 = c("-",
    "-"), V347 = c("-", "-"), V348 = c("-", "-"), V349 = c("-",
    "-"), V350 = c("-", "-"), V351 = c("-", "-"), V352 = c("-",
    "-"), V353 = c("-", "-"), V354 = c("-", "-"), V355 = c("-",
    "-"), V356 = c("-", "-"), V357 = c("-", "-"), V358 = c("-",
    "-"), V359 = c("-", "-"), V360 = c("-", "-"), V361 = c("-",
    "-"), V362 = c("-", "-"), V363 = c("-", "-"), V364 = c("-",
    "-"), V365 = c("-", "-"), V366 = c("-", "-"), V367 = c("-",
    "-"), V368 = c("-", "-"), V369 = c("-", "-"), V370 = c("-",
    "-"), V371 = c("-", "-"), V372 = c("-", "-"), V373 = c("E",
    "-"), V374 = c("H", "-"), V375 = c("D", "-"), V376 = c("I",
    "-"), V377 = c("G", "K"), V378 = c("T", "N"), V379 = c("C",
    "V"), V380 = c("Y", "L"), V381 = c("G", "E"), V382 = c("-",
    "A"), V383 = c("-", "S"), V384 = c("-", "F"), V385 = c("-",
    "D"), V386 = c("-", "S"), V387 = c("-", "E"), V388 = c("-",
    "F"), V389 = c("-", "K"), V390 = c("-", "E"), V391 = c("-",
    "I"), V392 = c("-", "N"), V393 = c("-", "R"), V394 = c("-",
    "V"), V395 = c("-", "L"), V396 = c("-", "F"), V397 = c("-",
    "G"), V398 = c("-", "E"), V399 = c("-", "E"), V400 = c("-",
    "G"), V401 = c("-", "Q"), V402 = c("-", "Q"), V403 = c("-",
    "Q"), V404 = c("-", "Q"), V405 = c("-", "G"), V406 = c("D",
    "E"), V407 = c("A", "E"), V408 = c("-", "S"), V409 = c("-",
    "Q"), V410 = c("-", "Q"), V411 = c("-", "E"), V412 = c("-",
    "G"), V413 = c("-", "V"), V414 = c("-", "I"), V415 = c("-",
    "V"), V416 = c("-", "E"), V417 = c("-", "L"), V418 = c("-",
    "E"), V419 = c("-", "R"), V420 = c("-", "E"), V421 = c("-",
    "Q"), V422 = c("-", "I"), V423 = c("-", "R"), V424 = c("-",
    "E"), V425 = c("-", "L"), V426 = c("-", "I"), V427 = c("-",
    "K"), V428 = c("-", "H"), V429 = c("-", "A"), V430 = c("-",
    "-"), V431 = c("-", "-"), V432 = c("-", "-"), V433 = c("-",
    "-"), V434 = c("-", "-"), V435 = c("-", "-"), V436 = c("-",
    "-"), V437 = c("-", "-"), V438 = c("-", "-"), V439 = c("-",
    "-"), V440 = c("-", "-"), V441 = c("-", "-"), V442 = c("-",
    "-"), V443 = c("-", "-"), V444 = c("-", "-"), V445 = c("-",
    "-"), V446 = c("-", "-"), V447 = c("-", "-"), V448 = c("-",
    "-"), V449 = c("-", "-"), V450 = c("-", "-"), V451 = c("-",
    "-"), V452 = c("-", "-"), V453 = c("-", "-"), V454 = c("-",
    "-"), V455 = c("-", "-"), V456 = c("-", "-"), V457 = c("-",
    "-"), V458 = c("-", "-"), V459 = c("-", "-"), V460 = c("-",
    "-"), V461 = c("-", "-"), V462 = c("-", "-"), V463 = c("-",
    "-"), V464 = c("-", "-"), V465 = c("-", "-"), V466 = c("-",
    "-"), V467 = c("-", "-"), V468 = c("-", "-"), V469 = c("-",
    "-"), V470 = c("-", "-"), V471 = c("P", "K"), V472 = c("C",
    "S"), V473 = c("-", "S"), V474 = c("-", "S"), V475 = c("-",
    "R"), V476 = c("-", "R"), V477 = c("-", "S"), V478 = c("I",
    "L"), V479 = c("S", "S"), V480 = c("-", "-"), V481 = c("-",
    "-"), V482 = c("-", "-"), V483 = c("-", "-"), V484 = c("-",
    "-"), V485 = c("-", "-"), V486 = c("-", "-"), V487 = c("-",
    "-"), V488 = c("-", "-"), V489 = c("-", "-"), V490 = c("-",
    "-"), V491 = c("-", "-"), V492 = c("-", "-"), V493 = c("-",
    "-"), V494 = c("-", "-"), V495 = c("-", "-"), V496 = c("-",
    "-"), V497 = c("-", "-"), V498 = c("-", "-"), V499 = c("-",
    "-"), V500 = c("-", "-"), V501 = c("-", "-"), V502 = c("-",
    "-"), V503 = c("-", "-"), V504 = c("G", "S"), V505 = c("P",
    "Q"), V506 = c("D", "-"), V507 = c("A", "-"), V508 = c("W",
    "-"), V509 = c("N", "-"), V510 = c("H", "-"), V511 = c("S",
    "D"), V512 = c("G", "E"), V513 = c("F", "P"), V514 = c("Q",
    "F"), V515 = c("C", "N"), V516 = c("-", "L"), V517 = c("-",
    "R"), V518 = c("-", "N"), V519 = c("-", "R"), V520 = c("Q",
    "K"), V521 = c("-", "P"), V522 = c("-", "-"), V523 = c("-",
    "-"), V524 = c("-", "-"), V525 = c("-", "-"), V526 = c("-",
    "-"), V527 = c("-", "-"), V528 = c("-", "-"), V529 = c("-",
    "-"), V530 = c("V", "I"), V531 = c("C", "Y"), V532 = c("-",
    "-"), V533 = c("-", "-"), V534 = c("-", "S"), V535 = c("-",
    "N"), V536 = c("-", "-"), V537 = c("-", "-"), V538 = c("-",
    "-"), V539 = c("-", "-"), V540 = c("-", "-"), V541 = c("-",
    "-"), V542 = c("-", "-"), V543 = c("-", "-"), V544 = c("-",
    "-"), V545 = c("-", "-"), V546 = c("-", "-"), V547 = c("-",
    "-"), V548 = c("-", "-"), V549 = c("-", "-"), V550 = c("-",
    "-"), V551 = c("-", "-"), V552 = c("-", "-"), V553 = c("G",
    "K"), V554 = c("L", "F"), V555 = c("V", "-"), V556 = c("K",
    "G"), V557 = c("S", "R"), V558 = c("W", "W"), V559 = c("N",
    "Y"), V560 = c("L", "E"), V561 = c("-", "I"), V562 = c("-",
    "T"), V563 = c("-", "P"), V564 = c("-", "E"), V565 = c("-",
    "K"), V566 = c("-", "N"), V567 = c("-", "P"), V568 = c("-",
    "-"), V569 = c("-", "-"), V570 = c("-", "-"), V571 = c("-",
    "-"), V572 = c("-", "-"), V573 = c("-", "-"), V574 = c("-",
    "-"), V575 = c("-", "-"), V576 = c("-", "-"), V577 = c("-",
    "-"), V578 = c("-", "-"), V579 = c("-", "-"), V580 = c("V",
    "Q"), V581 = c("G", "-"), V582 = c("D", "-"), V583 = c("F",
    "-"), V584 = c("V", "-"), V585 = c("V", "-"), V586 = c("Y",
    "-"), V587 = c("R", "-"), V588 = c("Q", "-"), V589 = c("E",
    "-"), V590 = c("-", "-"), V591 = c("-", "-"), V592 = c("-",
    "-"), V593 = c("-", "-"), V594 = c("-", "-"), V595 = c("-",
    "-"), V596 = c("-", "-"), V597 = c("-", "-"), V598 = c("-",
    "-"), V599 = c("-", "-"), V600 = c("-", "-")), class = "data.frame", row.names = c("UniRef90_A0A1V4QIJ3",
"UniRef90_A0A0L9V5F8"))

```

# Answer
> 3 votes
**Update**

It seems that part of the desired behavior is to write the sequences in the alignment *separately*. I am not sure that I understand why it is necessary to write them separately. Generally it is inadvisable to write many small files rather than one big file, based purely on issues of performance; you are likely to hit the I/O bound on performance very quickly, and that will be much slower.

One solution that has been suggested previously for this kind of problem is to write the whole alignment to a single file, and then post-process that filed in a second step to extract the sequences of interest. This could be as trivial as a `samtools faidx` call from the shell, or you could work with it in seqinr as suggested in the previous answer.

Of course, it seems like it would be simpler for most purposes to just have the big alignment file and use e.g. samtools indexing / faidx / biopython / whatever to extract the sequences of interest. But it seems that that is specifically not desired in this case.

**Original answer**

I would strongly suggest using the seqinr package, which is optimized for solving exactly these problems. See e.g. this similar question.

You will be able to read in an alignment (or whatever), cast it to as.matrix.alignment, do whatever your analysis is, then cast it back to as.alignment, which can then be passed to write.fasta to write the files.

Whenever possible it is preferable to work with such libraries, which handle these utility functions well, fast, and easy.

---
Tags: r, fasta, dataframe
---


# Question
Title: Can't locate object method "new" via package Bio::Tools::Primer3
1. Hi there I am still a novice at this. I'm trying to run this perl script called FullSSR https://sourceforge.net/projects/fullssr, by following the instructions in their 'README.TXT' https://sourceforge.net/projects/fullssr/files/.
2. I did everything as instructed by the app developer(s) which included among others, replacing the Primer3.pm file provided in the \*.tar.gz file into /usr/share/perl5/Bio/Tools/Run/; opening the config.txt and specifying the path of primer3\_core (file) and primer3\_config/ (folder); placing the sequence file I wanted to analyse into the FullSSR folder and executing the perl script using the command line: `perl 'FullSSR\(ver-1.1\).p; -i sequence.fasta -m 5`.
3. However, it was to no avail. Instead, I got the message "Can't locate object method "new" via package Bio::Tools::Primer3" (perhaps you forgot to load "Bio::Tools::Primer3?) at /usr/share/perl5/Bio/Tools/Run/Primer3.pm line 497.
4. I am going to try to install/load Bio::Tools::Primer3 with the cpan module (cpanm) and will see if that does it.
5. Perhaps someone out there could suggest a helpful (easy to understand) manual on CPAN?!

# Answer
I managed to figure out the answer to my question. I installed Bio::Tools::Primer3 by typing the following into terminal: '`sudo cpanm Bio::Tools::Primer3 --force`' to install it. It only works when I used 'sudo' and 'force'
> 1 votes

---
Tags: perl
---


# Question
Title: Clarification regarding (deepTools) bamCoverage output
Im using bamCoverage on a file and ive set the bin size to 1 but also normalizing using RPKM for defined regions using a bedgraph file. The output from bamCoverage however is in binned segments each with a (RPKM) score. Am I correct in interpreting then that

1 - the output bedgraph tells me the (total) RPKM score for that particular bin?

2 - to get the total RPKM score for the region defined in the bedgraph I need to sum all the RPKM scores given for that particular region?

# Answer
> 3 votes
1. Yes
2. The total RPKM for the region in the bedgraph file is the value already present in the bedgraph file. If you want to estimate the RPKM for a larger region composed of multiple bins then you can average the per-bin RPKM values. Note that adjacent bins with the same value are merged by deepTools, so you will need to account for this.

Note that RPKM is NOT a robust normalization and that I should remove it from deepTools so people stop using it.

---
Tags: rna-seq, deeptools
---


# Question
Title: Open chromatin region : nucleosome , mono-nucleosome , di-nucleosome
I m following this tutorial for my ATAC seq analysis. Here in of the section they show profile of open chromatin region plots.

```
**Plotting insert sizes with Greenleaf open, mono- and di-nucleosome profiles
This looks very similar to the image from the Greenleaf paper.

We can now annotate our nucleosome free (< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) as in the Greenleaf study.**

```

Now Since i have done the differential accessibility analysis. So from that information how do i find out region which are "nucleosome free (\< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) ".

My objective

* I would like to see this "nucleosome free (\< 100bp), mono-nucleosome (180bp-247bp) and di-nucleosome (315-437) " from my differential accessible region which was done from readcount and then I ran deseq2.
* My goal is not to plot the way they have done rather find those regions mentioned above and see it between different conditions.

Now based on the tutorial of that section where they do for each bam file.

Do i need to run the analysis on each bam file or each condition or Can I find out the above from my deseq2 analysis which has regions and those regions I have annotated it .

# Answer
> 1 votes
It's impossible to find nucleosome free and such regions without a BAM file. Your only options are:

1. Filter the BAM files to contain only the size-range of interest and rerun the counting and statistics.
2. Find or write a counting program that can track fragment sizes.

BTW, those size ranges could use some adjustments, you're throwing away all fragments between 100 and 180 bases using their ranges.

---
Tags: atac-seq
---


title: What's the meaning of Txn factor?

text:
# Question
Title: What's the meaning of Txn factor?
The transcription factor track in UCSC is named as Txn factor. Why Txn == Transcription?

# Answer
> 0 votes
Any other abbreviation is likely to work for "translation" as well and the x makes it sound a bit like transcription (txn would be pronounced something like "tixon").

---
Tags: ucsc, transcription-factors, transcriptome-regulation
---


# Question
Title: What is the best way to address the question of doublets and multiplets in a single cell RNA seq data set?
I have attached a histogram plot of the number of genes per cell in a single cell RNA seq data set of lung endothelial cells. I do not find a bimodal or multimodal distribution of the number of genes per cell within any particular condition. Doublets/multiplets should have an aberrantly high number of gene expression and if many, can lead to a bi or multimodal histogram plot of the number of genes detected per cell. Can I conclude that I have an insignificant presence of doublets or multiplets in my data set from this plot and hence continue with downstream analyses?

I am reluctant to use automatic doublet-finding algorithms because they have a high possibility to remove cells in a transitional differentiating phase. I also think that the default threshold of 2500 genes per cell given by Seurat might sometimes lead to the loss of good cells depending on the cell type.

Apart from trajectory analyses, are there any other analyses on single cell RNA seq data that are sensitve to doublets/multiplets?

Thanks in advance for your help.

# Answer
> 2 votes
Indeed there's no satisfactory doublet/multiplet identification methods as of 2020.

The most reliable way to remove doublets is still manual removal of doublet subclusters guided by marker genes.

Not only trajectory analysis can be impacted by doublets, but clustering and embedding can also be affected for obvious reasons.

---
Tags: rna-seq, scrnaseq, seurat, single-cell
---


# Question
Title: Can anyone point me towards a HiC Dataset with domains of high confidence?
is there a dataset online containing a HiC contact matrix and file with TAD boundary locations, which were generated, either by hand or automated, with high accuracy? Although I already used programs such as SpectralTAD or GMAP to generate TADs myself, a dataset which was computed with higher accuracy would be nice to check my project against.

Preferred attributes would be:

* human genome
* whole genome
* 10 kb resolution or higher

Thank you very much in advance!

edit: One dataset I found is http://chromosome.sdsc.edu/mouse/hi-c/download.html, which includes multiple 10 kb datasets with TADs called.

# Answer
> 1 votes
I guess your best chance is one of the in-situ HiC papers from the Lieberman-Aiden lab. There are not too many HiC datasets with 10kb-or-better resolution out there for human. I would check this one:

A 3D Map of the Human Genome at Kilobase Resolution Reveals Principles of Chromatin Looping

---
Tags: hi-c, chromosome-capture, domains
---


# Question
Title: How can I get or create a reference genome for Bacteria?
I am a computer engineer and nowadays trying to grasp some concepts of Bioinformatics particularly, reference genomes and genomic variants. My aim is to find the effect of sequence features on variant types that occur in bacteria genomes. I chose Escherichia coli as the bacteria for my study since there seem to be more publication and data about it and its genomic sequences.

When I made a research about variants, I found out that we need to have a reference genome for the comparison with other alleles in order to find variants. And, I even read that we can use genomic sequence of one bacteria as the reference but it won't represent the cohort in an accurate manner. I got the genomic sequence data for several strains of Escherichia coli bacteria that been sequenced before. However, most of these sequences are not complete ones. I filtered out the complete DNA genomes. When I checked the reports on these genomes, I learned that these are complete genomes but not reference genomes. Only the ones which have `RefSeq category: Reference Genome` in their report are reference genomes for Escherichia coli.

Learning that most of the data I got for different strains of Escherichia coli were not reference genomes, I tried to seek to find reference genomes and I got this website of NCBI ref sequences for prokaryotes. This website shows only two reference sequences which match with what I have find in all sequencing data I got before. So, I want to ask two questions

**First:** How can I get a reference sequence for other strains of E. coli?

**Second:** If I cannot get a reference sequence for other strains, can I use complete genome data as the reference one? When I check the complete genome data, I cannot understand their order because there are other data for same strains that show loci of sequences but in the complete genome, I cannot see any order of the loci. Even, I cannot see data of some loci in the complete genome which blows my mind more.

I know it may not be appropriate to ask two questions in a single post, however, these two questions seem to support each other. I will appreciate any of your help.

# Answer
> How can I get a reference sequence for other strains of E. coli?

You can use any sequence you want as reference for a variation analysis. It just depends on the question you want to answer.

So to keep it simple, your Question is:

**What are the variant region from a fresh sequenced *E. coli* K12 sample**?

To know the differences from your known strain *E. coli* K12. A reason for this could be that you work with the strain K12 through out several years. During your studies you changed the conditions for the bacteria, like adding antibiotics or change the growth media. So some of the *E. coli* would adapt to the antibiotic and become resistant to it or they adapt better to the new media. This happens because some nucleotide changed. With each new experiment you used the best adapted colonies from the previouse one.

So you sequence your bacteria and compare it with the *E. coli* K12 strain. In this case this is your reference. If you use another strain, like SE15 or SE11, as reference, your result would not give you any answer for your question.

If your used data and publications gives you the information which strain were used, then use the given strain as reference.
> 1 votes

# Answer
in genomics, organisms have a single reference genome which is updated( or not) over a period of time. E.coli has different strains and this means there might be reference genomes for each of the strain. If there is no reference genome available the its likely that in constructing the sequence of those strains they might have used any of those two as their reference. This information you can find by reading the paper that published those strains. Therefore I will advice you pick any of those two available reference genome as your reference and work with and make sure you indicate which one you used in your report or document.
> 0 votes

---
Tags: genome, assembly, variant-calling, reference-genome, bacteria
---


title: Tools for comparing/visualizing FASTAs?

text:
# Question
Title: Tools for comparing/visualizing FASTAs?
I have two FASTAs/assemblies of the same species, but the bases are somewhat different. I would like to explore this.

What tools/methods exist to compare two FASTAs and see the difference in structure between them? Are there visualization tools available, something like IGV?

I think users would need to do a multiple-sequence alignment between these FASTAs, and then check for the differences between them. However, I don't know if a "standard approach" exists to do this.

# Answer
You can try Mauve multiple genome aligner.
> 2 votes

# Answer
Depending on how long the sequences are, you could use the LAST Web Service:

http://lastweb.cbrc.jp/

It has an HTML view that shows base-level differences with colour-coding.

LAST \[download / command-line version\] can also be used for producing whole-genome comparison dot-plots.

Another thing I do for comparing assemblies is to chop one assembly into smaller pieces (e.g. 10kb), map the pieces to other other assembly to generate BAM output, then look at the BAM pileup file in a BAM viewer (e.g. Tablet).
> 2 votes

# Answer
For whole-genome align: Maybe you could use mummer? Mummer is good to align the whole genome of species, and it also has so many other tools to modify the result. (It can show the SNP.)

For short sequence: You could try mega X
> 1 votes

# Answer
Before you compare your genome assembly fasta, I would recommend you first do a reference guided scaffolding using the tool RagTag.

Once you are done you can compare the newly generated scaffolds (with other genomes or not) using the software BRIG. I made a video of how to use this software here

In addition you can also do a genome annotation on your scaffold using prokka. I made a video for that as well here

Good luck
> 1 votes

---
Tags: sequence-alignment, fasta, assembly, visualization
---


# Question
Title: How to deal with heterozygosity during polishing of genome assembly based on long reads?
All the long-read sequencing platforms are based on single-molecule sequencing which causes higher per-base error rates. For this reason a polishing step was added to genome assembly pipelines - mapping raw reads back to assembly and correcting details of the assembly.

I have decent PacBio RSII dataset of single individual genome of heavily heterozygous non-model species. Assembly went well, but when I tried to polish the assembly using quiver it could not converge over a couple of iterations and I bet it is because of too great divergence of haplotypes.

Is there any other way to polish a genome with such properties? For instance, is there a way to separate long reads by haplotype, so I could polish using one haplotype only?

# Answer
> 4 votes
A few possibilities:

**Falcon**

Try falcon and falcon-unzip. These are designed exactly for your problem and your data: https://github.com/PacificBiosciences/FALCON

**Not Falcon**

If you think you have assembled haplotypes (which seems reasonable to expect given enough coverage), you should be able to see the two haplotypes by just doing all pairwise alignments of your contigs. Haplotypes should show up as pairs of contigs that are MUCH more similar (even with a lot of between-haplotype divergence) than other pairs. Once you have all such pairs, you can simply select one of each pair to polish.

# Answer
> 3 votes
You could also have a go at Canu. It's designed for long-read assembly (both PacBio and Nanopore), although not specifically for complex population sequencing. It tries to strip a genome down into its unique components, and generates paths from those components that are well-supported from the reads.

With regards to polishing, it seems to be the case that polishing doesn't converge, and there will be a lot of variants that just oscillate between two possibilities. For me and at least one other person at London Calling this year, there was basically no gain in accuracy for polishing past the third iteration. I used my own error correction algorithm, but they used the more "standard" polishing with Pilon. For what it's worth, the nanopore WGS consortium used Racon for polishing their Canu assemblies.

# Answer
> 0 votes
There is a new tool called Hapo-G for diploid-aware polishing. Unfortunately, it seems it's designed for short reads only. It seems to be a more practical choice with respect to inferring the phase information on individual reads.

---
Tags: ngs, assembly, long-reads
---


# Question
Title: OpenAPI-Generator Python Client Instructions Examples
Trying to learn a little bit of Python and Biopython I ended up figuring out what a Rest API is when the RCSB PDB repository released its new Rest API.

I tried to retrieve info from the PDB repository using python3 requests,

then I managed to use OpenAPI-generator to generate what they call a 'Python3 Client'.

Now I am at loss I have this ApiClient class:

```
print dir():

 ['NATIVE_TYPES_MAPPING', 'PRIMITIVE_TYPES', '_ApiClient__call_api',
  '_ApiClient__deserialize', '_ApiClient__deserialize_date',
  '_ApiClient__deserialize_datetime', '_ApiClient__deserialize_file',
  '_ApiClient__deserialize_model', '_ApiClient__deserialize_object',
  '_ApiClient__deserialize_primitive', '__class__', '__delattr__',
  '__dict__', '__dir__', '__doc__', '__enter__', '__eq__', '__exit__',
  '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__',
  '__init__', '__init_subclass__', '__le__', '__lt__', '__module__',
  '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__',
  '__setattr__', '__sizeof__', '__str__', '__subclasshook__',
  '__weakref__', '_pool', 'call_api', 'client_side_validation',
  'close', 'configuration', 'cookie', 'default_headers',
  'deserialize', 'files_parameters', 'parameters_to_tuples',
  'pool', 'pool_threads', 'request', 'rest_client',
  'sanitize_for_serialization', 'select_header_accept',
  'select_header_content_type', 'set_default_header',
  'update_params_for_auth', 'user_agent']

```

and don't know how to use it. Googling for examples or searching for info on the OpenAPI-generator website doesn't seem to work for me, I find plenty of tutorials and examples on how to generate the client but none about how to use it.

The generated client class should help to better interact with the server as I understood (I could be wrong here). I am asking for helps and hints here because the OpenAPI webpage under Community points to stackoverflow under \[openapi-generator\] tag. Any tutorial, easy documentation that could help me understand the output of the generator?

# Answer
> 0 votes
OK fine, found useful information on READme.md generated by OpenAPI-generator itself

and tried `pdoc` to get html document of my OpenAPI-generator generated client Module

Moreover on stackoverflow How to generate API Documentation using openapi-generator? you can find how OpenApi Generator can do generation of html documentation from your yaml or json definition of your API.

For openAPI json api specification I found RapiPDF too : Web Component to generate PDF from OpenAPI Spec (json yml file).

---
Tags: python, public-databases, pdb, api
---


# Question
Title: biopython-motifs , How can I create motif from sequences of different length?
I want to create motif of different length kmers. Following code works for kmers of same length but not different length.What to do ?

```
from Bio import motifs

sequences = SeqIO.parse("myfile.fasta","fasta")
    instances = []
    for data in sequences:
      instance = data.seq
      instances.append(instance)

motif = motifs.create(instances)

```

# Answer
Motifs are not designed to be variable length in general since PSSMs/PVMs don't support insertions/deletions (indels). Instead, I suggest profile HMMs, such as those offered by HMMER3, which do support indels. You can parse HMMER3 output with Biopython via `Bio.SearchIO.HmmerIO`
> 3 votes

# Answer
All sequences represented in a Biopython motif object must be the same length. For example, a motif has a `__len__` property, which has to consist of a single value.

There are other methods that *do* support variable-length motifs, I suggest that you look into them, but unfortunately biopython seems to be limited to fixed-length motifs.
> 1 votes

---
Tags: python, motifs
---


# Question
Title: How to detect a mutation and predict its consequence?
I read lecture notes about mutations, and wondered what kind of algorithms are there to detect mutations?

How do we know if the gene is mutated or whether it is a sequencing error?

I saw this thread which is related, but I am not sure how to use the CIGAR information, e.g., is it 100% accurate? What is the underline mechanism to detect mutation? Is there a way to predict what the mutation will cause?

# Answer
> 10 votes
I'll follow up to the great answer from Kamil S Jaron:

Regarding predicting what the variant ("mutation" is a very loaded term) will do, there are a variety of tools. Chief among these are annovar and VEP. The general idea behind these is to classify the variants according to their overlap with genes, which codons they change (if any), how big that change is (e.g., changes in charge are more likely detrimental) and so on. One could also consider conservation, since if a position is highly conserved then changes in it are more likely to be detrimental.

If you really want to predict how a variant will change a protein's function then that usually requires prior knowledge about the proteins in question. Eventually someone will use machine learning to cull the literature and provide good predictions, but I haven't seen that yet.

# Answer
> 9 votes
**Part 1 : how to detect mutations**

The keywords you are searching for are "variant calling". Basically you have to map sequencing reads to a reference genome (or gene) and then estimate for each position of the genome if the observed difference of mapped reads and the reference is more likely a sequencing error or a mutation (in genomic glossary - variant).

Popular tools for variant calling are GATK, FreeBayes or bcftools (previously part as Samtools package).

The question you linked asks for quick alternatives to simply looking for variants. Indeed, you can just visualise the mapped reads to the reference sequence and see if the variant is there or not. CIGAR is just a notation of read alignment used in sam files (files with mapped reads), you can find good explanation of CIGAR strings here.

The follow-up about "How to estimate an effect of mutation" is in @DevonRyan 's awesome answer.

# Answer
> 3 votes
## Part 3

> I thought I'd add a **part 3** . Having found via VEP a truncation (missense or frameshift) or a missense, how do you find what are the consequences at the protein level?

The effect a mutation has on a protein can be diverse and how much imbalance in function can the cell's equilibrium absorb depends on the component. As a result you get recessive, dominant and embryonically lethal mutations. In this figure are shown some of the effects, divided into core and surface mutations:

## Core variants

Core mutations are the most black-and-white to analyse, but require a protein structure available. A severe truncation and a highly deleterious missense will result in a protein that will likely be degraded. The difference in folding Gibbs free energy (∆∆G, where $\Delta \Delta G = \Delta G\_{\rm{mutant}} – \Delta G\_{\rm{WT}}$)) can be calculated for a missense, but how large this must be to result in a certain percentage being degraded depends on factors that cannot be calculated analytically. Alternatively, the Gibbs free energy of ligand binding is affected resulting in loss of catalysis.

To calculate the ∆∆G effect (stability), there are many tools:

* Online SDM uses ML to estimate the ∆∆G quickly (do note that positive kcal/mol normally is a loss of stability, but here the sign is inverted so positive is a gain in stability
* Rosetta toolkit. A highly flexible panel of movers (samplers) can calculate ∆∆G
* FoldX
* MD runs
* Missense3D assesses a variant based on several structural criteria

In the case of ligand binding modelling as above but using the ligand bound structure and calculating the binding energy as the difference between bound and unbound (far far away). For altered catalysis, modelling as above but using the transition state and assessing the difference in the transition hill can be done. However, an active site variant is very likely to affect either ligand binding, transition state binding (catalysis) or product binding that it's surely deleterious.

This is all (mostly) straightforward. If the domain is a regulatory domain however, there is loss of inhibition, which results in a gain of function (dominant variant).

## Surface variants

Then stuff gets messier when you get to surface variants and these can cause dominant effects. However, if you are close to a phosphorylation, ubiquitination or other PTM site, as identified in phosphosite plus, you likely have a case of deregulation. For linear motifs, ELM is useful, but requires a lot of literature scanning to confirm any find. Alternative a loss of an interface may break the protein function, however, what binds where can only be ascertained by available structures, cross-linking mass spec or in some cases well defined Western experiments, so this data is rare. For this type of variant there are two pipelines:

* Miscast collates information on the protein
* VENUS (beta) calculates ∆∆G and finds any relevant close-by PTM or gnomAD variants \[disclaimer: this is mine\].

## Footnote

As the variants can have a variety of effects, none of these *score* (unlike the genetics parts), but instead try to *inform* as to what may be going on.

---
Tags: variant-calling, gene
---


title: ATAC seq : Plotting insert sizes

text:
# Question
Title: ATAC seq : Plotting insert sizes
My tutorial link

```
library(GenomicAlignments)

atacReads <- readGAlignmentPairs("MACS2_peak_call/ALL_ATAC_data/HSC1.bam", param = ScanBamParam(mapqFilter = 1,
            flag = scanBamFlag(isPaired = TRUE, isProperPair = TRUE), what = c("qname",
        "mapq", "isize"), which = GRanges("chr20", IRanges(1, 63025520))))
# length(atacReads)
atacReads

atacReads_read1 <- GenomicAlignments::first(atacReads)
insertSizes <- abs(elementMetadata(atacReads_read1)$isize)
head(insertSizes)

library(magrittr)
library(dplyr)
library(ggplot2)
fragLenPlot <- table(insertSizes) %>% data.frame %>% dplyr::rename(InsertSize = insertSizes,
    Count = Freq) %>% mutate(InsertSize = as.numeric(as.vector(InsertSize)),
 Count = as.numeric(as.vector(Count))) %>% ggplot(aes(x = InsertSize, y = Count)) +
  geom_line()

fragLenPlot + theme_bw()

fragLenPlot + scale_y_continuous(trans = "log2") + theme_bw()

```

How to i specify the Grange to span the total chromosome. Now I use only "chr20" as shown in the example.

# Answer
> 1 votes
You can exclude the `which` statement and get the entirety of all of the chromosomes. Note however that R is a very bad platform to use for this, since you end up reading everything into memory (I hope the BAM file isn't huge).

Note that you're seriously reinventing the wheel. There are a number of packages that will directly produce the plot you want, such as `bamPEFragmentSize` from deepTools.

---
Tags: atac-seq
---


# Question
Title: Nvidia Parabrick fq2bam pipeline error - No such file or directory
I am trying to implement Nvidia-parabrick fq2bam pipeline to process my WES dataset. The reason for opting this pipeline is the huge number of samples (1004) which take many months to process for somatic mutation identification. THe GPU accelerated pipeline will be more helpful to get somatic mutations in a quick manner. After pulling parabrick docker image from here, the following command works fine and return the options :

> docker run -it nvcr.io/hpc/parabricks:v2.5.0 -h

But whenever I give the following command to run fq2bam pipeline, I am getting an error of `No such file or directory`:

> docker run -it nvcr.io/hpc/parabricks:v2.5.0 fq2bam --ref ucsc.hg19.fasta --in-fq sample\_1.fastq.gz sample\_2.fastq.gz --knownSites 1000G\_phase1.indels.hg19.sites.vcf.gz --out-bam sample\_dups\_gpu.bam out-recal-file sample\_recal\_gpu.txt --tmp-dir /tmp

The error is as follows:

```
Traceback (most recent call last):
File "/parabricks/run_pipeline.py", line 7, in <module>
sys.exit(PB.pb_main())
File "PB.pyx", line 846, in PB.pb_main
File "/parabricks/pbargs.py", line 1287, in getArgs
return PBRun(sys.argv)
File "/parabricks/pbargs.py", line 392, in __init__
self.runArgs = getattr(self, args.command)(argList)
File "/parabricks/pbargs.py", line 447, in fq2bam
self.addToParser(bsmb_parser_sysgroup, sysOptionGenerator().allOptions)
File "/parabricks/pbargs.py", line 50, in __init__
PBOption(category="sysOption", name="--num-gpus", default=GetNumGPUs(), typeName=int, helpStr="Number of GPUs to use for a run")
File "/parabricks/pbutils.py", line 48, in GetNumGPUs
output = subprocess.check_output(["nvidia-smi", "-L"], universal_newlines=True)
File "/usr/lib/python2.7/subprocess.py", line 567, in check_output
process = Popen(stdout=PIPE, *popenargs, **kwargs)
File "/usr/lib/python2.7/subprocess.py", line 711, in __init__
errread, errwrite)
File "/usr/lib/python2.7/subprocess.py", line 1343, in _execute_child
raise child_exception
OSError: [Errno 2] No such file or directory

```

As I am using the docker image as it is,it is not clear to me which file is missing. Kindly suggest. Thanks.

# Answer
My query was answered by NVIDIA Parabrick developer team. The installation instruction are available at https://forums.developer.nvidia.com/uploads/short-url/t5uXsZgdMXx0G2wDcYamFl9sO2x.pdf.
> 1 votes

---
Tags: sam, genomics, fastq
---