# Question
Title: Finding matched genes for each genomic range
I have called structural variants and now I have

```
Sample  Type    Chrom1   Start1        End1 Strand1
t_005   DEL 1   28677074    28677079    +
t_005   INS 1   186383731   186383738   +
t_005   DEL 1   245902589   245902590   +
t_005   DEL 2   56345643    56345645    +
t_005   DEL 3   59766214    59766217    +
t_005   DEL 3   60270545    60270548    +

```

I have an annotation file like

```
> head(sv[,1:5])
  ENSG00000233440 13 23708313 23708703   HMGA1P6
1 ENSG00000207157 13 23726725 23726825    RNY3P4
2 ENSG00000229483 13 23743974 23744736 LINC00362
3 ENSG00000229483 13 23743974 23744736 LINC00362
4 ENSG00000252952 13 23791571 23791673  RNU6-58P
5 ENSG00000235205 13 23817659 23821323  TATDN2P3
6 ENSG00000232849 13 93708910 93710179 LINC00363
>

```

How I could relate each gene name from annotation file to my genomic range file?

> EDITED

```
Sample  Type    Chrom1  Start1  End1    Strand1 Chrom2  Start2  End2    Strand2
t_005   DEL 1   28677074    28677079    +   1   28706324    28706329    -
t_005   INS 1   186383731   186383738   +   1   186383731   186383732   -
t_005   DEL 1   245902589   245902590   +   1   246007384   246007385   -
t_005   DEL 2   56345643    56345645    +   2   56345829    56345830    -
t_005   DEL 3   59766214    59766217    +   3   59785325    59785328    -
t_005   DEL 3   60270545    60270548    +   3   60705591    60705594    -
t_005   DUP 3   60744344    60744345    -   3   60806306    60806307    +
t_005   DUP 3   60882109    60882112    -   3   60934155    60934158    +
t_005   DEL 3   77973091    77973093    +   3   77991398    77991400    -
t_005   DEL 3   95972796    95972800    +   3   95974511    95974515    -

```

# Answer
> 3 votes
You received an answer from Kamil for a similar question here: https://bioinformatics.stackexchange.com/a/8532/650

You can adapt the code they gave you to compare `Chrom1 == sv[2] & Start1 >= sv[3] & End1<=sv[4]` to find events that happen within a gene. For events spanning genes, you're going to need to write a collapsing function for the gene names and a different way to compare overlaps instead of `Start >= sv[3] & End <= sv[4]`.

---
Tags: r, structural-variation, wgs
---


title: Getting this information for GISTIC2

text:
# Question
Title: Getting this information for GISTIC2
Long time I am struggling how to provide ***GISTIC2*** required input like `Num_Probes` and `Segment_Mean`

I have called somatic copy number by scatngs and I have copy number output file (attached);

https://www.dropbox.com/s/r6nef2v1j2di6q6/my.copynumber.caveman.csv?dl=0

This is few lines of scatngs copy number output

```
1,1,13116,665401,2,1,3,1
2,1,668374,96324627,2,1,2,1
3,1,96326095,96419709,2,1,1,0

```

In scatngs reference paper says in this file we have

```
Segment number

Chromosome

Start position (origin-1)

End position (origin-1)

Major copy number—normal

Minor copy number—normal

Major copy number—tumor

Minor copy number—tumor

```

scatngs also returns another output like below

https://www.dropbox.com/s/7i28vooixdp06c2/WTSI-OESO\_121\_1pre.copynumber.txt?dl=0

```
Chromosome  Position    Log R   segmented LogR  BAF segmented BAF   Copy number Minor allele    Raw copy number
rs62635286  1   13116   -0.0290741822638953 0.245246748556562   0.25    NA  3   1   2.35806652809599
rs75454623  1   14930   0.402869048014209   0.245246748556562   0.3095  NA  3   1   4.14031506947036

```

I am wondering how I can get `Num_Probes` and `Segment_Mean` From this files?

People say

```
Segment_Mean = log2(tumour copy number / 2)

```

So, does that mean I should add up `Major copy number—tumor` and

`Minor copy number—tumor` to get tumour copy number? for getting Num\_Probes what should I do?

***EDITED***

> people says
>
> scatngs uses read counts at SNPs to generate calls and so the locations of the SNPs is what I need to generate the number of markers in each segment. (eg each SNP = 1 marker).
>
> Do they mean, for instance if the position of SNP rs62635286 is 13116, so the num\_prob is 13116?

Any help please?

Thanks a lot in advance

# Answer
Depends on the software to generate the segmentation file;

For example I have used scatngs for segmentation that for sample returns

```
Segment number
Chromosome
Start position (origin-1)
End position (origin-1)
Major copy number—normal
Minor copy number—normal
Major copy number—tumor
Minor copy number—tumor

```

By log2(Major copy number—tumor + Minor copy number—tumor /ploidy of each sample) you would get segment mean

For number of probes

```
The number of probes depends on where the data came from.
1) From SNP arrays:  Number of array probes within each segment.
2) From exon capture data:   Probes are the ends of exons and so the number of probes are equal to the number of exons  X 2.
3) From whole genome sequencing:  Probes are the boundaries of bins used to generate the copy number data and so the number of probes are equal to the number of bins +

```

For example in my case for scatngs

```
Each SNP is a marker, and so you have to count how many of these SNPs are on each segment

```

I have asked this question in

Extracting the number SNP in each range
> 0 votes

---
Tags: wgs, copy-number
---


# Question
Title: GSEA - Gene Set Enrichment Analysis - in R studio
I have been trying to conduct a classic GSEA using the broad institute software and can not figure out how to compile it in R for the life of me. I am aware that there exists other packages in R that are *almost* the same, such as `fgsea` however, I would like to use this one first and then compare it with these other methods.

I have already generated a ranked list for my contrast of interest and save it as a `.rnk` file. Here is what I tried to do:

```
gseaCallH <-
  paste("java -Xmx2G -cp C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out C:/Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7",sep="")

testGSEA <- system(command = gseaCallH, intern = TRUE, ignore.stderr = TRUE)

```

When I run the `system` command I get the following error: `had status 1`

At first I thought maybe it had something to do with JAVA not being properly installed however, I have installed it and still nothing. Maybe the installation was not done correctly? I have tried running it directly from cmd (I'm running windows 10) and got the following error:

```
Error occurred during initialization of VM
Could not reserve enough space for 2097152KB object heap

```

Can anyone advise on how to proceed?

# Answer
> 3 votes
I have answered my own question however, it is specific to windows users.. If you get the same error from the `cmd` as I did, you must manually put in environment variables related to Java. Do as the user Mohamed did in this SO post: https://stackoverflow.com/questions/9303889/error-occurred-during-initialization-of-vm-could-not-reserve-enough-space-for

Following this, you **must avoid** using paths such as `C:/absolute_path` and instead use `/absolute_path`, if you do the former, you will get the exception: `java.net.MalformedURLException: unknown protocol: c`

In the end, after manually entering the environment variables and running the following command: `java -Xmx2G -cp /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea2-2.2.4.jar xtools.gsea.GseaPreranked -param_file /Users/omk/Desktop/BIOINFORMATIC/MsigDB/GSEAPrerankedParameters.txt -rnk /Users/omk/Desktop/BIOINFORMATIC/MsigDB/gsea_treatFU_min.rnk -gmx /Users/omk/Desktop/BIOINFORMATIC/MsigDB/h.all.v6.2.symbols.gmt -rpt_label gsea_treatFU_min.rnk -out /Users/omk/Desktop/BIOINFORMATIC/MsigDB/C7` it has worked.

Where the command can be executed directly from the command line or in R using the `system` function.

**Note**: A '#' separated by a *space* must be present at the beginning of your `.rnk` file for the program to work. Therefore, `.rnk` files must have the following format:

`# Genes Coeffecients
TSPAN6 0.2345676` etc. for every gene.

# Answer
> 2 votes
Are you calling java with sufficient memory? Try increasing `-Xmx2G` to `-Xmx8G`, provided your PC has 8 GB RAM.

---
Tags: r, rna-seq, gsea, java
---


# Question
Title: Alignment for predicting DNA hybridization?
I am currently working on a Computer Science project where we are trying to build a large set of orthogonal single-stranded DNA sequences. The goal would be to ensure that when put in solution, the ssDNA do not interact with each other to form double-stranded DNA.

I hope that this makes sense, but my knowledge biological understanding is developing because I have a CS background. Please excuse any inadvertent errors.

As one of the first steps, we used NUPACK to predict the percent of the resulting AB complex, given two single-stranded sequences A, B in concentration 1μM each. We would use the concentration of AB to classify if the strands are likely to form double-stranded DNA.

However, we cannot use NUPACK for large sets, as we would need to do the above check against each of the other sequences, which is a quadratic problem. For a set of size 30,000 we would have to do 450 million NUPACK computations which is simply unfeasible.

So I turned to the idea of approximating NUPACK. In this paper, it is claimed that

> Ideally, the thermodynamics of hybridization between a candidate probe and potential off-target sites would be modeled in silico and employed as a means of identifying probe oligos likely to only bind their intended targets in a given set of reaction conditions. Although powerful utilities such as NUPACK (52– 54) are capable of performing such simulations, the limited throughput of these programs renders a direct thermodynamic approach impractical for genome-scale probe design. However, we hypothesized that features in rapidly calculated data such as alignment scores may be predictive of thermodynamic behavior and could therefore serve as a proxy for the information that would be produced by thermodynamic simulations.

The alignment scores are then used to predict the interaction. bowtie2 is used in the paper for alignment purposes. However, the context is different as their computation is done on genome information. Very sensitive settings of bowtie2 seem to be slow (hours) on my laptop and result in enormous files (tens of GB). I used less sensitive settings as a compromise but it seems to negatively affect the LDA.

I am looking at alternative applications for alignment, but there are many alternatives (some of the more promising ones so far, in my opinion: SNAP, MMseqs2, Genoogle). Would any of these, or others not mentioned, be better suited for my purpose?

Also any other advice or suggestions are welcome!

# Answer
As a rough and ready and very easy to perform test I would put the entire universe of dsDNA through Primer3 (details below). This could be done by "splicing" your probe of interest onto each "chunk" of dsDNA you pass through primer3 and restrict the primer search to the 20bp primer you are interested in. In other words you are forcing the program to consider the internal and trans-hybridisation of every possible dsDNA against your probe. If you probe is not recovered from a given search you would flag this on a subsequent parsing script and then that "chunk" of DNA would be assessed further. All the components required are in the link below. You would need to provide your own wrapper, but its fairly easy.

Primer3 will assess the internal thermodynamics of your probe and it has loads of options to assess the full range of behaviour, for example 3' hybridisation (which could be extended because you want 5' to 3'). It would be easy to construct a parameter grid around all the parameters of interest. Ultimately it would provide an additional layer of certainty, and everyone knows and trusts this algorithm (its also written in C, so its fast).

The description of all the components is here , my specific purpose for customising the pipeline is different.

---

Bowtie2 will take a lot of RAM, I am not sure I would use it for this purpose, particularly if RAM is limiting (its used for genome assemblies against a reference). I'd ask the users whether Blat would be useful, because that is friendlier to the calculation you are performing. Anyway there some heavy duty Bowtie2 users around and their input would be useful.
> 1 votes

---
Tags: sequence-alignment, fasta, blast, phylogenetics, thermodynamics
---


# Question
Title: WGCNA error codes in Network Construction & moduel detection
I am using WGCNA and code has worked with no problems with some of my datasets however with most recent attempt when completing gene network & identification of modules. I get following error code...

```
net__30= blockwiseModules(t5, power=6,
                              TOMType = "signed", minModuleSize=30,
                              reassignThreshold =0, mergeCutHeight = 0.25,
                              numericLabels = TRUE, pamRespectsDendro = FALSE,
                              saveTOMs = TRUE,
                              saveTOMFileBase = "PBCreponseTOMB",
                              verbose= 3)

```

I come up with error code saying:

```
            Number of items to replace is not a multiple of replacement length.

```

Does anyone please know what this indicts to? I have looked on all the WGCNA tutorial pages but cannot find this error message anywhere.

Thank you kindly

# Answer
Looks like a bug in the code which I will try to hunt down and fix. In the meantime, I would suggest playing with maxBlockSize argument to blockwiseModules. Try increasing it as much as your available RAM allows (see the paragraph "A second word of caution concerning block size" on page 6 of WGCNA tutorial I, section 2c, for some guidance on how to set maxBlockSize). Hope this helps.
> 2 votes

---
Tags: r, wgcna
---


# Question
Title: How to get Uniprot and swissprot cross references using Biopython
I am struggling to use biopython to gather cross references such as GO annotation from Uniprot/swissprot text files ("`DR`" in the .txt files) using biopython.

From the source code:

```
def _read_dr(record, value):
    cols = value.rstrip(".").split('; ')
    record.cross_references.append(tuple(cols))

```

This leads me to think that the following code should work, as this is how other features like annotation and features work.

```
for record in SeqIO.parse(filename, "swiss"):
    print(record.cross_references) #does not work, should return a list of DR rows
    print(record.features) #works, returns a list of parsed FT rows

```

# Answer
> 3 votes
I believe the attribute you should be looking for is `dbxrefs`:

```
[x for x in record.dbxrefs if x.startswith('GO:')]

```

Apparently `SeqIO` handles uniprot text files differently than the `SwissProt` module:

```
from Bio import SwissProt
record = SwissProt.read(open('single_record_from_uniprot.txt'))
record.cross_references

```

So you should change you code to use the `SwissProt` module to get the behavior you expect

---
Tags: python, uniprot
---


# Question
Title: Demographic model for admixed African Americans
I am fairly new to genomics and population genetics. And I am reading this article 'Demographic history and rare allele sharing among human populations':

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3142009/pdf/pnas.201019276.pdf

In this article they propose a demographic model in Figure 4 on page 11986. However, I am confused how does African Americans fit in this demographic model that are known to have both African and European ancestry. And if I am to consider African American population than what is the split time between the two populations, time they merge into admixed population and the admixed proportion. My apologies if it's a really basic question but, any insights would be appreciated.

# Answer
Just to further my answer, the problem with molecular dating within a species (I don't do humans \[excepting immunology\], but it's the same theory), is it requires an assumption of clonality (bifurcating trees). Mendelian species recombine towards panmixia, which will make any attempt to date via de novo point mutation, appear far younger than they really are, or for panmixia, incalculable (it breaks the central assumption). Dating through admixture requires a lot of extra data and/or assumptions for example the generation time, putative admixture per generation, and shifts in population dynamics. Maybe you can calculate some of those requirements but alot will remain assumptions.

I have not read the paper you suggested, because I'm not involved in human genetics \[outside immunology\], but it is possible the authors may (and I stress may) have simplified an assumption that could be difficult to estimate. It is not clear what the objective of the paper was.

Quantitative analysis is extremely important in modern population genetics/ evolutionary theory because historically it was overlooked and resulted in lots of hypothese which were difficult to test.
> 1 votes

---
Tags: phylogenetics
---


title: How is BLAST's nr database created?

text:
# Question
Title: How is BLAST's nr database created?
Is there a paper or web page describing the procedure for creating the nr database used by NCBI's BLAST implementation?

I presume it's some type of clustering, but I'm curious about how exactly sequences are condensed into non-redundant representatives.

# Answer
Did a little more searching and found the answer in the README on BLAST's ftp site: ftp://ftp.ncbi.nlm.nih.gov/blast/db/README

```
6. Non-redundant defline syntax

The non-redundant databases are nr, nt and pataa. Identical sequences are
merged into one entry in these databases. To be merged two sequences must
have identical lengths and every residue at every position must be the
same.  The FASTA deflines for the different entries that belong to one
record are separated by control-A characters invisible to most
programs. In the example below both entries Q57293.1 and AAB05030.1
have the same sequence, in every respect:

>Q57293.1 RecName: Full=Fe(3+) ions import ATP-binding protein FbpC ^AAAB05030.1 afuC
[Actinobacillus pleuropneumoniae] ^AAAB17216.1 afuC [Actinobacillus pleuropneumoniae]
MNNDFLVLKNITKSFGKATVIDNLDLVIKRGTMVTLLGPSGCGKTTVLRLVAGLENPTSGQIFIDGEDVTKSSIQNRDIC
IVFQSYALFPHMSIGDNVGYGLRMQGVSNEERKQRVKEALELVDLAGFADRFVDQISGGQQQRVALARALVLKPKVLILD
EPLSNLDANLRRSMREKIRELQQRLGITSLYVTHDQTEAFAVSDEVIVMNKGTIMQKARQKIFIYDRILYSLRNFMGEST
ICDGNLNQGTVSIGDYRFPLHNAADFSVADGACLVGVRPEAIRLTATGETSQRCQIKSAVYMGNHWEIVANWNGKDVLIN
ANPDQFDPDATKAFIHFTEQGIFLLNKE

Individual sequences are now identifed simply by their accession.version.

For databases whose entries are not from official NCBI sequence databases,
such as Trace database, the gnl| convention is used. For custom databases,
this convention should be followed and the id for each sequence must be
unique, if one would like to take the advantage of indexed database,
which enables specific sequence retrieval using blastdbcmd program included
in the blast executable package.  One should refer to documents
distributed in the standalone BLAST package for more details.

```

Landed on that README from this question on biostars.org: https://www.biostars.org/p/217456/

**Edit**

In that same README file is some information on the origin of the sequences in the non-redundant sets:

```
+-----------------------+-----------------------------------------------------+
|File Name              | Content Description                                 |
+-----------------------+-----------------------------------------------------+
nr.gz*                  | non-redundant protein sequence database with entries
                           from GenPept, Swissprot, PIR, PDF, PDB, and RefSeq
nt.gz*                  | nucleotide sequence database, with entries from all
                          traditional divisions of GenBank, EMBL, and DDBJ;
                          excluding bulk divisions (gss, sts, pat, est, htg)
                          and wgs entries. Partially non-redundant.

```
> 2 votes

# Answer
The Refseq team and also the NCBI resource coordinators team publish a new paper every few years, so check out the many papers (e.g. here or here), but to answer your 2nd question, non-redundancy here is (I think) defined very strictly as proteins that are identical in terms of sequence and length, so the clustering is trivial, without the need for a sophisticated clustering algorithm as required to detect more remote homologs.
> 4 votes

---
Tags: blast, makeblastdb
---


# Question
Title: Database for germline copy number variations in cancer
I am interested in looking at germline copy number variations in individuals that are at high risk of developing cancer. Are there any databases where I can look, if a CNV seen in our test case has been seen/reported earlier in a patient with cancer? Thanks.

# Answer
> 2 votes
You can find CNV (and other structural variant) information in:

* DECIPHER

> DECIPHER is used by the clinical community to share and compare phenotypic and genotypic data. The DECIPHER database contains data from 29603 patients who have given consent for broad data-sharing; DECIPHER also supports more limited sharing via consortia. Have a look at the numbers.
>
> Anyone can browse publicly-available patient data on DECIPHER and request to be put in contact with the responsible clinician. Why? Because sharing benefits everyone.
* ExAC

> The Exome Aggregation Consortium (ExAC) is a coalition of investigators seeking to aggregate and harmonize exome sequencing data from a wide variety of large-scale sequencing projects, and to make summary data available for the wider scientific community.
>
> The data set provided on this website spans 60,706 unrelated individuals sequenced as part of various disease-specific and population genetic studies. The ExAC Principal Investigators and groups that have contributed data to the current release are listed here.

ExAC has included CNV calls since August 2016.
* ClinVar

> ClinVar is a freely accessible, public archive of reports of the relationships among human variations and phenotypes, with supporting evidence. ClinVar thus facilitates access to and communication about the relationships asserted between human variation and observed health status, and the history of that interpretation. ClinVar processes submissions reporting variants found in patient samples, assertions made regarding their clinical significance, information about the submitter, and other supporting data.
* Database of Genomic Variants

> The objective of the Database of Genomic Variants is to provide a comprehensive summary of structural variation in the human genome. We define structural variation as genomic alterations that involve segments of DNA that are larger than 50bp. The content of the database is only representing structural variation identified in healthy control samples.
>
> The Database of Genomic Variants provides a useful catalog of control data for studies aiming to correlate genomic variation with phenotypic data. The database is continuously updated with new data from peer reviewed research studies. We always welcome suggestions and comments regarding the database from the research community.

Information from all of these is also displayed in the Structural Variant browser of VarSome*. While you can't search by patient status, so you can't use VarSome to know the clinical history of any sample with the CNV, you can at least look up a specific gene and see the structural variants listed for it. For example, for TP53, you would see:

---

* Disclaimer: I work for the company behind VarSome.

---
Tags: variant-calling, variation, cancer, cnv
---


# Question
Title: Sort vcf by contig and position within contig
Due to tabix constraints, I need to sort a vcf so that contigs and then positions within contigs occur in numerical order in the vcf. I don't know if the following will sort positions within contigs, but I tried this:

```
cat in.vcf | vcf-sort -c > out.vcf

```

but was informed:

```
Old version of sort command installed, please run without the -c option.

```

I've read that can cause problems here. Further, my contigs are't numeric, but are something like this: 'Erumb1\_s00000030'.

I see two possible solutions. (1) I am on Mac OS Mojave, and if there were a way to install the a new unix-style sort, I could use it. I tried homebrew and didn't find anything. (2) there are likely some bash commands that would do it. If anyone has an idea in either of these directions, I would appreciate it.

# Answer
You could simply use bcftools sort for it:

```
$ bcftools sort input.vcf > output.vcf

```

If you really want to use bash only, you can do this:

```
$ grep "^#" input.vcf > output.vcf
$ grep -v "^#" input.vcf| sort -k1,1V -k2,2g >> output.vcf

```

The first command will write the header information to the new vcf file. The second will sort by contig name and position and append the result to the new vcf file.

fin swimmer
> 5 votes

# Answer
Use a subshell. First print the header and then the variant lines, sorting by chromosome and position:

```
(grep ^"#" input.vcf; grep -v ^"#" input.vcf | sort -k1,1 -k2,2n) > sorted.vcf

```
> 2 votes

---
Tags: vcf, bash, vcftools
---


# Question
Title: Merging transcriptomes coming from different experiments
I'm planning to build a transcriptome by pooling all existing transcriptomes in SRA for a non-model species (which has no reference genome) to study differentially expressed genes and the like. It might worth mentioning that the transcriptomes I want to pool mostly come from different tissues and sampling sites. These transcriptomes would come from five SRA BioProjects plus one project from my laboratoryfor a total of 102 single transcriptomes coming from about 20 individuals (plus three pooled larvae transcriptomes). I'm especially interested in differentially expressed genes in different tissues, on genes specific to different tissues and on the UTRs of differentially expressed genes.

My idea is to build a common reference transcriptome against which I'd re-assemble the individual transcriptomes, and only then I'd get differential expression data etc.

I'm looking for literature (or personal experience) dealing with assemblies of transcriptomes / comparing transcriptomes coming from different datasets in order to understand whether this is a **good idea or not** and get a grasp of what **technical and biological problems** I could face. I'm (vaguely) aware of the batch effect and of one way to remove it, but I'd like to read papers dealing with my same problem (different SRA databases, without reference genomes). Up until now I found just projects dealing with microarray data.

# Answer
> 2 votes
As mentioned in the comments, your biological question is going to be vitally important in determining the best way to tackle this issue. However, there are some things that you should generally keep in mind.

1. **Most *de novo* assemblers will perform kmer normalization.** This means that the merged RNA-Seq data will be normalized across all libraries. This is almost certainly not what you want because you'll be randomly selecting reads from all libraries without regard to their treatments. So, if you're going to allow for kmer normalization, you should normalize prior to merging the libraries so that you have equal representation of the reads by library. If you don't normalize, you're going to have a very fragmented assembly, and it will take a long time to generate.
2. ***De novo* assemblers using short reads assume all identified exons with junction support may be used to generate splice variants with equal probability.** I'm not entirely sure how to best phrase that, but here's an example. Suppose in a cold-treated library that Gene A generates a transcript using exons 1-2-3, but a heat-treated library generates transcripts using 3-4-5 and 2-3-5. It may be that merging the libraries results in transcripts generated containing various permutations of those exons (e.g. 1-2-3-5, 2-3-4, etc.) that aren't generated in any of the biological settings but can no longer be distinguished. This will lead to false positive transcripts which will detrimentally affect expression quantification.
3. **Merging assemblies is difficult.** So, to be true to each treatment, it makes more sense to generate individual assemblies (which will also result in library-specific kmer normalization). StringTie essentially does this using reference genomes. You assemble using each library (all bio-reps for a treatment), and then merge after the fact and re-quantify. This is the best approach. However, it can be difficult and computationally taxing. Some programs like CD-HIT (specifically cd-hit-est) are typically used in these cases. Keep in mind that CD-HIT will drop sequences with high similarity, not merge them. This may lead to loss of some uniquely assembled sequence. Other approaches involve using programs like CAP3 (a DNA assembler) to treat the transcripts like longer reads for assembly. Depending on the parameter choice, this might work, but you will risk generating chimeric transcripts.

How you move forward with the assembly process will again depend on your biological questions. I would tend to promote that you tackle the questions one at a time rather than using a single transcriptome to address them all, but there may be circumstances in which merging assemblies makes more sense.

In your specific case, I would perform multiple assemblies using a variety of the approaches described above and then annotate the assemblies (ex. using Trinotate) to see how their quality metrics compare. I would discourage generating a single transcriptome and then re-assembling individual transcriptomes from that assembly -- the single transcriptome may have misassemblies that would then propagate to the individual assemblies, and this process would generally be messy/difficult.

---
Tags: database, differential-expression, transcriptome, literature-search
---


title: Reading a txt in R with weird header

text:
# Question
Title: Reading a txt in R with weird header
I am trying to read a txt file in R but always I am getting

```
Error in scan(file = file, what = what, sep = sep, quote = quote, dec = dec,  :
  line 113 did not have 10 elements

```

I noticed that in the header of my file I have something like

```
 0000775 €   g¿´�€   4_®™00000016125 13404727020 017402  0                                                                                ustar   devons01                        stlab-icgc

```

Even by removing these characters I am not able to read this file

This is link of my file

https://www.dropbox.com/s/rpc5jgti5aa433b/txt.txt?dl=0

Could somebody please heap me in reading this file in R?

# Answer
> 2 votes
Delete the corrupt portion preceeding `CFS` in the header and it works with with `read.delim()`.

# Answer
> 3 votes
We can skip the first row:

```
read.delim("txt.txt", skip = 1, header = FALSE)

```

**Or** we can get the *"uncorrupted"* version of the table, it looks like it is from this paper, Table S4.

---
Tags: r
---


# Question
Title: Calculate the occurrence of motif in sequences (including overlaps) using Python
I need to calculate the occurrences of a motif (including overlaps) in sequences (motif is passed in the first line of standard input and FASTA sequences in subsequent lines). The sequence name starts with \>, and after whitespace is just a comment about the sequence that needs to be neglected. The input of program is like:

```
 AT
 >seq1 Comment......
 AGGTATA
 TGGCGCC
 >seq2 Comment.....
 GGCCGGCGC
 GCCGGCGGG

```

The output should be:

```
seq1: 2
seq2: 0

```

I decided to save the first line as a motif, strip the comment from sequence name, join lines of sequence in one line and save sequence names (keys) and sequences (values) in a dictionary. I also wrote a function for motif\_count and want to call it on dictionary values and then save it in a separate dictionary for final output. However, it does not work the way I want it as the method motif\_count does not produce any result when called on my d values.

Here is my code:

```
 #!/usr/bin/env python3

 import sys

 sequence = sys.stdin.readlines()
 d = {}
 temp_genename = None
 temp_sequence = None
 motif = sequence[0]

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
    if start <0:
        return count
    count += 1

if sequence[1][0] != '>':

    print("ERROR")

exit(1)

for line in sequence[1:]:
    if line[0] == '>':
       temp_genename = line.split(' ')[0].strip()
       tmp_genename = temp_genename.replace('>', '')
       temp_sequence = ""

else:

      temp_sequence += line.strip()

d[tmp_genename] = temp_sequence
for key in d:
    motif_count(motif, key)

```

# Answer
> 1 votes
I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):

```
#!/usr/bin/env python3

import sys

sequence = sys.stdin.readlines()
d = {}
temp_genename = None
temp_sequence = None
# The lines returned by readlines() contain end-of-line characters,
# this is why strip() is needed
motif = sequence[0].strip()

def motif_count(motif, temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1

if sequence[1][0] != '>':

    print("ERROR")

    exit(1)

for line in sequence[1:]:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            d[temp_genename] = temp_sequence
        temp_genename = line.split(' ')[0].strip()
        temp_genename = temp_genename.replace('>', '')
        temp_sequence = ""

    else:

        temp_sequence += line.strip()

# This should be done each time a new sequence is found,
# not only at the end
d[temp_genename] = temp_sequence
for key in d:
    # We need to use the result
    # motif_count(motif, key)
    # We want to count the motif in the sequence, not in its name
    count = motif_count(motif, d[key])
    print(f"{key}: {count}")

```

Testing it:

```
$ cat tmp.fa
AT
>seq1 Comment......
AGGTATA
TGGCGCC
>seq2 Comment.....
GGCCGGCGC
GCCGGCGGG
$ cat tmp.fa | ./count_motif.py
seq1: 2
seq2: 0

```

Here is another way to do, using `valmap` to make a dict of the motif counts from the dict of sequences, and with a few other programming style changes:

```
#!/usr/bin/env python3

import sys
from cytoolz import valmap

motif = sys.stdin.readline().strip()

# We can make this a function of only one parameter
# because motif is fixed
def motif_count(temp_sequence):
    count = 0
    start = -1
    while True:
        start = temp_sequence.find(motif, start +1)
        if start <0:
            return count
        count += 1

seq_dict = {}
temp_genename = None
temp_sequence = None

for line in sys.stdin:
    if line[0] == '>':
        # Record the previous one in the dictionary
        if temp_genename is not None:
            seq_dict[temp_genename] = temp_sequence
        temp_genename, *_ = line[1:].strip().split()
        temp_sequence = ""
    else:
        temp_sequence += line.strip()

msg = f"""File format might be wrong"
motif: {motif}
temp_genename: {temp_genename}
temp_genename: {temp_sequence}"""

assert temp_sequence is not None and len(temp_sequence), msg

# Record the last sequence
seq_dict[temp_genename] = temp_sequence

# Make a new dict where the values are the motif counts
# valmap needs the function motif_count to take only one parameter
counts_dict = valmap(motif_count, seq_dict)
# Print the results, iterating over dict items
for (name, count) in counts_dict.items():
    print(f"{name}: {count}")

```

# Answer
> 3 votes
This is a good task for Biopython:

```
from Bio import SeqIO

with open('input_file') as f:
    motif = next(f).rstrip().lower()
    for r in SeqIO.parse(f, 'fasta'):
        print(r.id, ':', r.seq.lower().count_overlap(motif))

```

Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the standard file format.

# Answer
> 0 votes
I like @Chris\_Rands tidy use of OO, objects are really important per se.

There appears a paradox with your assignment that you might want to refer back to your supervisor/tutors.

If "AT" is the motif, this cannot overlap with another "AT" motif, therefore the request for "overlapping motifs" makes this part of the code superfluous.

It would be better expressed if the motif was "ATA" for example. Thus if the sequence was ATATATA the motif is present 3 times, but only twice if the motif was contiguous.

More generally the code essentially searches for AT microsatellite motifs and formally performed using e.g. tandem repeat finder.

---
Tags: python, phylogenetics, homework
---


title: How about BioMoose?

text:
# Question
Title: How about BioMoose?
This is a silly post, but hope it makes a point. We have BioPerl (which never impressed me, in fact there isn't even a BioPerl tag), BioPython (which admittedly is impressive), and even now BioRust.

Why not BioMoose?

Moose is Perl5s answer to objects and is really impressive in its simplicity and ease of use. Before Moose Perl5 objects, were/are, how can I say, make evolutionary chaos theory look hip and relevant. Moose changed that completely, but stacking Moose code into BioPerl is never gonna work.

# Answer
> 2 votes
I fear this entire thread is going to veer dangerously close to opinionation. But regardless of how one feels about the Perl language and its various libraries and communities, I think we can objectively say that the era of Perl's dominance in the bioinformatics community has passed. A lot of tools and libraries implemented in Perl are still widely used, and there is a tremendous amount of legacy Perl code around. The ability to read and write Perl remains a valuable skill.

But in 2019, few people (with a choice) are learning Perl or implementing new tools and frameworks in Perl. Python fits snuggly into the niche once occupied by Perl, and its native OO support is (I dare to say) objectively superior to any *post hoc* OO support Perl 3rd-party libraries provide.

The broader web/tech community moved on from Perl to Python (and Ruby and others) even before the bioinformatics community did, and this broader community's offerings (IPython, Jupyter, scipy, numpy, pytest, etc.) make Python an even more attractive alternative to Perl for bioinformatics work.

In the meantime, the "hadleyverse" ecosystem of R packages for data manipulation make R much more capable of the types of data pre-processing tasks for which an R programmer might have previously used a language like Perl or Python.

In short, Perl is waning in popularity (if not relevance), and I would not expect much enthusiasm from the bioinformatics community for a new framework built around a relatively obscure Perl extension.

---
Tags: perl
---


# Question
Title: BUSCO does not find buscos in reference sequence of C. elegans (technical problem)
I am working on a GNU/Linux machine. I have installed BUSCO using conda. I have downloaded the nematode dataset `nematoda_odb9.tar.gz`here. I downloaded the reference sequence of *C. elegans* from Ensembl. I then ran busco using the following command:

`run_busco --in Caenorhabditis_elegans.WBcel235.dna.toplevel.fa --cpu 1 --out busco --lineage_path ./nematoda_odb9/ --mode genome --force > busco.o 2> busco.e`

tblastn seems to have run properly. `run_busco/blast_output/tblastn_busco.tsv` contains lots of hits for many of the buscos.

However, I have 982 "Missing BUSCOs".

The file that contains logs and errors related to augustus (run\_busco/augustus\_output/augustus.log) contains:

```
Will create parameters for a EUKARYOTIC species!
creating directory ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/ ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_parameters.cfg ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_weightmatrix.txt ...
creating ...anaconda/envs/busco3_0_2/config/species/BUSCO_busco_1395218643/BUSCO_busco_1395218643_metapars.cfg ...
The necessary files for training BUSCO_busco_1395218643 have been created.
Now, either run etraining or optimize_parameters.pl with --species=BUSCO_busco_1395218643.
etraining quickly estimates the parameters from a file with training genes.
optimize_augustus.pl alternates running etraining and augustus to find optimal metaparameters.

```

It does not seem to contain any errors.

I have the following error:

`Failed to open sequence file run_busco/augustus_output/extracted_proteins/EOG091H07ZC.faa.1` In fact, the folder `run_busco/augustus_output/extracted_proteins/` is empty.

Anyone would know why I have no extracted proteins?

# Answer
Most likely your Augustus installation is not functioning properly. Did you install Augustus from source or using conda? I believe the conda version of Augustus is known to have some issues.

Something worth trying: Running a complete busco install in a VM. This image has all dependencies installed and properly configured.
> 2 votes

---
Tags: software-usage, hmmer, busco
---


title: Complexheatmap update issue

text:
# Question
Title: Complexheatmap update issue
After i updated my complexheatmap library i can;t use my old annotation to label sample this is the image how i used to get something like this

Code for the below image

```
require(ComplexHeatmap)
require(circlize)
require(cluster)
df <- read.csv("heatmap_folder/NORMAL_EPI_TF1.txt",header = TRUE,sep = "\t")

head(df)
tail(df)
dim(df)
heat <- t(scale(t(df[,3:ncol(df)])))
head(heat)

rownames(heat) <- df$Gene
Sample = c(rep("HSC",4),rep("CMP",4),rep("GMP",4),rep("Mono",4))

anno <- as.data.frame(Sample)
categories <- data.frame(anno)
colnames(anno) <- c("Condition")

ColAnn <- HeatmapAnnotation(df=anno, which="col", col=list(Condition=c("HSC"="black","CMP"="brown",
                                                                       "GMP"="maroon1","Mono" = "darkgreen")))
RowAnn <- data.frame(df$Family)
colnames(RowAnn) <- c("Gene family")
RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")

boxAnnCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                    pch=".", size=unit(2, "mm"), axis=TRUE, axis_side=NULL, axis_gp=gpar(fontsize=10)),
                               annotation_width=unit(c(1, 6.5), "cm"))

boxAnnRow <- rowAnnotation(boxplot=row_anno_boxplot(heat, border=TRUE, gp=gpar(fill="#CCCCCC"),
                                                    pch=".", size=unit(3, "cm"), axis=TRUE, axis_side="top", axis_gp=gpar(fontsize=12)),
                           annotation_width=unit(c(3), "cm"))

myCol <- colorRampPalette(c("navyblue", "white", "red"))(100)
myBreaks <- seq(-2, 2, length.out=100)
hmap <- Heatmap(heat,
                name="Z-score",
                col=colorRamp2(myBreaks, myCol),
                heatmap_legend_param=list(color_bar="continuous",
                                          legend_direction="vertical", legend_width=unit(5,"cm"),
                                          title_position="topcenter", title_gp=gpar(fontsize=10, fontface="bold")),

                #Split heatmap rows by gene family
                split=df$Family,
                #split = NULL,

                #Row annotation configurations
                cluster_rows=TRUE,
                show_row_dend=FALSE,
                row_dend_side = "left",
                #row_title="Transcript", #overridden by 'split' it seems
                row_title_side="left",
                row_title_gp=gpar(fontsize=8,fontface="bold"),
                row_names_gp = gpar(fontsize = 8),
                show_row_names=TRUE,
                row_names_side="right",
                row_title_rot=0,

                #Column annotation configuratiions
                cluster_columns=FALSE,
                show_column_dend=TRUE,
                column_title="Samples",
                column_title_side="top",
                column_title_gp=gpar(fontsize=10, fontface="bold"),
                column_title_rot=0,
                column_names_gp = gpar(fontsize = 20, fontface="bold"),
                show_column_names=TRUE,

                #Dendrogram configurations: columns
                clustering_distance_columns="euclidean",
                clustering_method_columns="complete",
                column_dend_height=unit(10,"mm"),

                #Dendrogram configurations: rows
                clustering_distance_rows="euclidean",
                clustering_method_rows="complete",
                row_dend_width=unit(15,"mm"),

                top_annotation_height=unit(0.5,"cm"),
                top_annotation=ColAnn,

                bottom_annotation_height=unit(2, "cm"),
                bottom_annotation=boxAnnCol)
hmap + RowAnn

```

The code remains the except the data

Now not sure may be im messing up the annotation so that i can;t see sample label and the family annotation

Since this was the error

> "Error: `axis_gp` is removed from the arguments. Use `axis_param =
>  list(gp = ...)` instead." which i modified a bit

Then i modified my code to this

```
RowAnn <- data.frame(df$Family)
colnames(RowAnn) <- c("Gene family")
RowAnn <- HeatmapAnnotation(df=RowAnn, col=colours, which="row")

ColAnn <- HeatmapAnnotation(df=anno, which="col", col=list(Condition=c("HSC"="black","CMP"="brown",
                                                                       "GMP"="maroon1","Mono" = "darkgreen")))

boxplotCol <- HeatmapAnnotation(boxplot=anno_boxplot(heat,
  border=TRUE, gp=gpar(fill="#CCCCCC"), pch=".",
  size=unit(2, "mm"), axis=TRUE, axis_param = list(side = "left") ,
  annotation_width=unit(c(1, 5.0), "cm"), which="col")
)

hmap <- Heatmap(heat, name = "expression",col=colorRamp2(myBreaks, myCol),
                split=df$Family,
                #split = NULL,

                show_row_names = TRUE, show_column_names = TRUE, cluster_rows = TRUE,
                cluster_columns = TRUE, show_column_dend = TRUE, show_row_dend = TRUE,
                row_dend_reorder = TRUE, column_dend_reorder = TRUE, clustering_method_rows = "ward.D2",
                clustering_method_columns = "ward.D2", width = unit(100, "mm"),

                #top_annotation_height=unit(1.0,"cm"), top_annotation=colAnn,

                bottom_annotation=boxplotCol)

draw(hmap+RowAnn, heatmap_legend_side="left", annotation_legend_side="right")

```

my fig

I would be glad to know what im doing wrong as i can;t she condition getting annotated

my data

```
dput(df)
structure(list(Family = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("EPI",
"TF"), class = "factor"), Gene = structure(c(31L, 23L, 32L, 40L,
26L, 34L, 14L, 17L, 39L, 21L, 15L, 30L, 9L, 1L, 33L, 2L, 5L,
11L, 18L, 16L, 7L, 27L, 8L, 19L, 22L, 12L, 24L, 36L, 6L, 20L,
28L, 41L, 43L, 44L, 45L, 47L, 48L, 49L, 50L, 13L, 17L, 38L, 42L,
21L, 25L, 35L, 37L, 10L, 11L, 4L, 46L, 3L, 29L), .Label = c("BABAM1",
"BAP1", "BRF2", "CASZ1", "CBX8", "CCDC17", "CEBPA", "CUX2", "DAPK3",
"DEAF1", "DOT1L", "EBF1", "ERF", "HCFC1", "HDAC7", "HES6", "HMG20B",
"INO80E", "IRF8", "KLF16", "MBD3", "MEF2D", "MTA3", "MXD4", "PIN1",
"PIWIL4", "PREB", "REPIN1", "RFX2", "SETD1A", "SFMBT1", "SFMBT2",
"SMARCA4", "SMYD2", "SNAPC4", "SREBF2", "TCF20", "TCF7", "TFPT",
"UBE2A", "ZBTB17", "ZFHX2", "ZNF142", "ZNF319", "ZNF414", "ZNF446",
"ZNF469", "ZNF581", "ZNF646", "ZNF76"), class = "factor"), HSC1 = c(1.38256100990102,
6.03508433561645, 8.53718694596963, 9.7579, 4.9089643439748,
6.08652168925513, 10.85362, 5.791612, 4.717658, 5.034092, 7.6573073927514,
8.214206, 6.087519, 7.795752, 9.147491, 9.16065391580272, 5.475244,
6.839048, 4.307532, 0.4008883181761, 3.681525, 5.39357937244362,
0.048702874829245, 4.884674, 9.42704096616912, 0.331449426359954,
8.70131868583511, 9.860997, 5.96926894038996, 5.122815, 9.84364182958174,
4.99646644487682, 8.64613509396529, 6.221252, 4.788327, 4.221885,
6.904238, 7.027222, 6.90290115677703, 5.910114, 5.791612, 7.85232611449929,
4.30633834185685, 5.034092, 6.587824, 3.77111448331312, 4.816565,
5.83621984608683, 6.839048, 5.38002542796833, 2.499083, 4.07296733289583,
7.384605), HSC2 = c(3.3250440950498, 5.69710569315867, 7.53535864842004,
9.46294, 4.63015366085456, 6.27268560573257, 10.61808, 5.936903,
5.754081, 5.17955, 7.34546746135814, 8.156995, 4.959736, 7.1724,
9.519366, 8.399964626033, 3.79417, 6.580323, 3.006501, 0.465543296427136,
3.854056, 6.19801925234426, 1.35661238539194, 2.520579, 10.1179535946398,
0.438637606305892, 8.02299949380521, 9.599955, 5.14592215971399,
5.362818, 9.62067506274524, 5.81177950264385, 8.6388082408386,
6.580456, 5.507812, 4.647564, 5.933962, 7.036506, 6.61438422233323,
7.081847, 5.936903, 7.25061985280206, 6.14119015507351, 5.17955,
5.51992, 4.7980903060972, 4.688224, 5.75480261400296, 6.580323,
4.61480421865767, 2.473632, 2.58901029457056, 7.141953), HSC3 = c(4.5609308176212,
6.31335736623603, 8.08532370258223, 9.954938, 5.33971874563167,
6.31519305036765, 8.686445, 3.652774, 1.705097, 3.590782, 5.82598582988908,
7.160368, 2.730684, 6.544555, 8.308113, 7.28529068473834, 2.181592,
6.153712, 3.769699, 0.686412239668329, 3.181293, 5.81660138585872,
0.513314099277583, 3.807037, 9.40253894125251, 0.789263521770297,
6.25643885795079, 8.4667, 3.59885151800005, 1.959458, 8.16946518123644,
3.78892554507724, 8.17931637873862, 6.871673, 2.207662, 2.586419,
4.611411, 6.342816, 6.24601538582422, 4.504617, 3.652774, 6.1448035804027,
4.32474896196878, 3.590782, 4.283432, 3.97377694799887, 5.081613,
4.67353303760197, 6.153712, 2.95856351440761, 1.797669, 2.61144324906345,
6.927889), HSC4 = c(4.33336593087016, 5.51218211320577, 7.98927439884575,
10.02818, 4.72376095157983, 7.06408657043682, 10.6532, 6.332262,
4.356266, 5.063746, 7.40605558817077, 7.943625, 4.577773, 7.114057,
10.10959, 8.39439253099313, 0.843862, 7.275198, 3.125764, 0.461081947661647,
4.914624, 6.47822557945205, 0.151047922603591, 4.233378, 9.76778876120242,
1.52153581933298, 7.74186897540428, 9.848883, 4.27034862337376,
5.461201, 9.20573611283668, 5.57160015049504, 8.90499660789483,
5.972766, 4.37926, 4.914465, 6.288035, 6.667292, 6.48347655570048,
6.525983, 6.332262, 6.67910931337338, 5.93921207229708, 5.063746,
5.391649, 5.14382963917955, 3.682687, 6.17864627473077, 7.275198,
4.55820389425563, 2.464959, 3.36552973589941, 7.155267), CMP1 = c(6.87345497292835,
5.51412758368083, 7.31352446967914, 10.55862, 6.25444074607741,
7.17818981714724, 11.09789, 6.353011, 4.176257, 5.988031, 8.45885235860671,
8.239023, 4.831419, 8.439534, 9.741489, 9.12973766558991, 1.63498,
7.616455, 4.40003, 3.35836219130892, 6.782445, 8.09320436383016,
0.034087182831375, 7.328077, 10.1348328626404, 3.77197918672699,
8.22105235434344, 9.586611, 3.96736165755976, 6.120908, 9.66405823977144,
5.5410623284247, 7.98439042490943, 5.617969, 4.389489, 3.320773,
7.259381, 6.342179, 7.80157963446663, 6.00359, 6.353011, 4.88520066048698,
4.57336386910792, 5.988031, 7.534274, 4.22034258843835, 4.407149,
6.04554630191099, 7.616455, 5.5286242078284, 2.938658, 6.05835585199437,
7.254343), CMP2 = c(4.57632794651802, 3.78653056181795, 7.38323687896953,
10.18887, 5.42618113653324, 7.14753074148422, 11.10573, 6.444083,
5.364945, 5.56755, 7.87006213457315, 8.183806, 5.547153, 7.408484,
10.41511, 8.5734105703595, 3.332623, 7.927959, 3.044397, 2.23975823913519,
5.948004, 6.8523382099618, -0.208208601046198, 2.502613, 9.61958525630466,
0.083082776643226, 7.96816766770755, 9.658819, 4.28029266262047,
6.064432, 9.51407062066367, 6.0002887227803, 8.84860518364821,
6.63929, 5.418921, 4.246354, 6.265858, 6.726515, 6.59033513636581,
6.599091, 6.444083, 6.42468503826279, 5.39157032731654, 5.56755,
5.975197, 5.49143188875569, 4.882971, 6.3269787223519, 7.927959,
4.52014312274917, 3.02202, 3.32107044636493, 7.234948), CMP3 = c(5.85154481845215,
4.70657054545281, 7.81269743162002, 10.51041, 6.02525567584962,
7.6061510389796, 8.938033, 4.185122, 3.453336, 4.250061, 6.95084523398581,
7.016165, 3.711631, 6.691113, 9.291314, 6.79870749458327, 2.556527,
6.640901, 2.472737, 0.713500223112626, 4.610029, 7.09501723097758,
0.554991742557483, 3.663828, 8.64957961123419, 0.830920666997163,
5.03120897448563, 8.193863, 3.66385830275847, 2.514419, 8.57920717695102,
3.68106350062846, 7.99067978896729, 4.917504, 1.705394, 2.940584,
4.935565, 4.679952, 5.45403468978382, 4.097788, 4.185122, 4.7017134210513,
3.36996680968093, 4.250061, 4.841303, 4.16528140537129, 4.311709,
4.34280333979994, 6.640901, 3.46696931369128, 1.836416, 1.2954688478789,
6.759544), CMP4 = c(5.70607301130076, 4.90995364796149, 8.07668148492955,
10.49127, 6.12272998130959, 8.1070578791571, 11.14272, 6.069435,
3.77716, 5.905144, 7.80284615350651, 8.193236, 5.266133, 7.994962,
10.62971, 8.76622729912917, 1.728943, 8.018743, 2.64787, 3.17362412630739,
6.675791, 8.17425862958478, 1.3062705649262, 8.305058, 9.60100662458438,
3.64083547966437, 7.03151917955059, 9.697494, 4.22387311586196,
5.872986, 9.48040032459999, 5.9054376735385, 9.0993222961645,
6.641411, 4.042785, 4.175351, 6.040673, 7.025595, 6.8701243794046,
6.474802, 6.069435, 5.98875791029547, 5.44536133412665, 5.905144,
6.253757, 5.55713735563886, 3.183721, 6.23569788279688, 8.018743,
4.14168543985438, 3.140228, 2.85944735485205, 7.19532), GMP1 = c(6.00003380096822,
6.45782274084263, 10.1386873750581, 11.18212, 4.99574155426759,
8.22548270576452, 11.80921, 6.139842, 4.716938, 6.208399, 7.25702470622416,
9.216353, 4.459075, 7.604192, 10.214, 9.71979154806924, 4.045693,
6.892321, 4.36733, 4.54209088244969, 4.363115, 7.42142574723769,
8.16658165526266, 13.13801, 11.7030896033115, 1.38480651276548,
7.70549837966537, 10.51655, 5.6799276696584, 7.117451, 10.4613436109924,
7.51060888343217, 9.45360701936786, 11.02082, 4.741664, 4.37359,
6.757818, 7.101331, 6.51838172412493, 6.234017, 6.139842, 7.16147421243415,
5.80674368005962, 6.208399, 6.860234, 4.69772329684801, 5.404108,
6.45675195710298, 6.892321, 5.30230900977469, 4.155601, 4.36874064027233,
7.42003), GMP2 = c(5.74747276487971, 5.97553824097351, 8.84397150139003,
10.54739, 7.01942050935297, 8.32418895763336, 11.28484, 6.100887,
5.371866, 5.782092, 6.81459497341517, 8.300889, 4.573727, 7.155416,
10.55908, 8.53584780642442, 3.288921, 7.981419, 3.343707, 1.78160205563733,
7.106303, 6.97534663269409, 4.30909405726115, 10.94426, 10.9546928815164,
6.74862372900741, 7.1870575783895, 9.650893, 4.16099729988147,
6.067064, 9.97780340347263, 5.92239739476255, 9.00840651100957,
5.952921, 5.034178, 4.255218, 5.533772, 6.671744, 6.58968157482356,
6.572507, 6.100887, 6.06061888187195, 4.67987580614607, 5.782092,
5.976481, 5.57605980178547, 4.822138, 5.9343704738438, 7.981419,
4.0036610231274, 2.709241, 2.35084397835661, 7.122418), GMP3 = c(6.86258223460792,
6.20078322375291, 9.61446552902509, 10.80572, 7.36372917472567,
8.3826076358515, 11.15008, 5.201591, 4.123945, 5.045783, 6.24154676477513,
8.276453, 4.263376, 7.519041, 10.51785, 8.30200623104098, 2.495567,
7.857207, 2.962492, 1.93915883789869, 7.249551, 7.84810925448229,
5.84602007858824, 11.77414, 11.1027524778454, 5.63020389504886,
6.50727637405781, 9.456897, 3.75134646474614, 5.932853, 10.1039559568546,
4.9400044296399, 9.31472753734898, 6.039605, 3.336343, 3.975937,
5.736387, 6.905666, 6.30344093479224, 5.426484, 5.201591, 5.47026133533826,
4.69976751703472, 5.045783, 5.749752, 5.23937658540781, 5.154187,
6.34478978856221, 7.857207, 3.50805939441741, 2.657933, 2.78283797334871,
7.091083), GMP4 = c(6.05880869638876, 6.28654396576068, 9.96479543345143,
10.73202, 7.213420739218, 8.81669595423024, 11.31576, 5.679182,
4.155314, 5.673872, 6.43309255133785, 8.167664, 4.513147, 7.768381,
10.68025, 8.7651559260891, 2.485242, 8.41621, 3.457858, 2.07304058721841,
6.412872, 8.09668562132081, 5.97986866984983, 12.08623, 11.28032878957,
6.57392060927003, 6.41585754599334, 9.749641, 3.66579392385564,
5.878271, 9.75282317403084, 5.77880844881496, 9.31419691100924,
5.048616, 3.784422, 4.082659, 6.015668, 6.789027, 6.44947035574528,
6.177861, 5.679182, 5.40028170642952, 4.40100191946731, 5.673872,
6.189491, 5.32217020480945, 3.106884, 5.98621706786421, 8.41621,
3.507010215991, 2.810708, 2.26597487448895, 6.822586), Mono1 = c(5.2563684147039,
6.45633325704225, 8.35683375157621, 10.4683, 6.89578871222261,
6.89453589831897, 10.76922, 5.79888, 3.938016, 4.625604, 8.14811443428627,
8.89819, 4.565085, 7.582637, 9.690787, 9.26841770880937, 2.680727,
7.256856, 5.506176, 1.73004263683502, 7.168846, 7.34133258574044,
2.04394933289888, 10.88391, 10.8831286687177, 2.76098817217713,
7.62084289021189, 9.376334, 5.89374731019185, 5.059797, 9.58406921227736,
6.10083260496288, 8.90601784817498, 8.926934, 4.514644, 4.671412,
6.824876, 7.160325, 7.12646156474527, 6.029051, 5.79888, 7.45127447114964,
6.58727075697179, 4.625604, 6.726376, 4.54492833429885, 5.482923,
5.75789405204644, 7.256856, 5.01661617666884, 3.912327, 4.21481739920369,
7.311136), Mono2 = c(5.85774218629914, 5.24626015596457, 7.87393184409552,
10.74763, 6.1115477644807, 8.59285304183197, 11.30387, 6.487999,
5.641296, 5.605635, 6.83338037544917, 8.158746, 5.031193, 7.791207,
10.12365, 8.83146107303257, 3.839879, 7.995034, 3.864015, 1.95600811707846,
6.547528, 7.15767104450585, 1.13087340949285, 10.60541, 10.6662569861529,
2.45950067788357, 7.80075639569934, 9.772267, 4.10768104209722,
6.039518, 9.45126490091168, 5.43166366925388, 8.53913541490463,
6.852944, 5.210281, 3.909477, 6.10825, 6.759114, 6.10127101835935,
6.394382, 6.487999, 6.77443082143061, 4.65095032762504, 5.605635,
6.292733, 5.06179443435785, 4.687693, 5.97873952456647, 7.995034,
4.53685791636727, 2.036068, 3.79175025763666, 7.80933), Mono3 = c(5.79480048210947,
5.93841021449181, 8.4456428811783, 10.96997, 5.90368127873194,
7.57633522462406, 10.77521, 5.357294, 3.631277, 5.573695, 6.02937219504517,
7.823033, 4.89389, 7.134801, 9.942904, 8.29297682562671, 2.376142,
7.942619, 2.938838, 2.58663314475541, 5.473615, 6.54736080917345,
1.50818476189011, 10.39009, 9.81896143227971, 2.61179572528994,
7.11294273571537, 9.321912, 4.59467076052555, 4.846162, 9.0889984655259,
5.5580711224014, 8.60684080318993, 6.807437, 3.781213, 3.349656,
5.352098, 6.732808, 6.18873397274794, 5.580769, 5.357294, 6.02896733254353,
4.09719798968178, 5.573695, 5.576249, 4.76264681239467, 4.599397,
5.85107919870538, 7.942619, 4.46991267319193, 2.359072, 2.49361466378781,
7.168415), Mono4 = c(5.79136484133828, 5.22885638262284, 8.66188302349651,
10.80097, 6.02816710441609, 8.28180453176157, 10.95579, 5.606963,
4.781397, 6.146859, 7.35941757150855, 8.279649, 5.035288, 8.131193,
10.46082, 8.7442177262891, 2.904981, 8.009644, 2.855958, 2.31248706248342,
7.171729, 7.85215015030227, 1.5674747074644, 10.34766, 10.4041048711624,
1.7198811173707, 7.7219124608883, 9.652986, 4.80881131147327,
5.631832, 9.37540653744023, 5.93520028734304, 8.83460813839155,
7.145424, 4.609511, 4.550807, 6.666205, 6.956604, 6.66113029874611,
6.54366, 5.606963, 5.47022342513485, 4.59138757680696, 6.146859,
5.969648, 5.20105550006015, 3.857023, 6.22204033269913, 8.009644,
4.44000233834592, 2.854773, 2.88863557103663, 6.788677)), class = "data.frame", row.names = c(NA,
-53L))

```

# Answer
I see that you created the `ColAnn` variable with the Condition annotation but did not include it in the `draw` call (or even in the `Heatmap()` function). Could that be the problem here?

EDIT: Oops, just noticed that this post is *really old*. Doesn't have an answer, so I guess I'll let mine stay.
> 1 votes

---
Tags: r, visualization, heatmap
---


# Question
Title: How to learn/get started bioinformatics?
I know this is a broad question and I'm sorry if this is not the right place to post (this is my first post ever), but is there some sort of roadmap for bioinformatics?

I'm an undergraduate student and I have a strong background in statistics, computer science, and molecular biology so I feel like I should be able to get around, but every time I talk or read about bioinformatics I feel lost. It's such a vast area that I don't even know where to begin or what even exists. Is there some kind of all-encompassing textbook or roadmap for me to follow?

Most textbooks I could find expend way to much time explaining really basic stuff like what is a cell and how to program and I feel they fail to convey the bigger picture of bioinformatics. Things like what exists in the area, what are the current technologies, conventions, and challenges that we face. What should I know in order to call myself a bioinformatician?

# Answer
> 2 votes
The formal answer is O'Reilly "Bioinformatics programming, using Python". Check out the classic O'Reilly texts on "Mastering Perl for bioinformatics", and possibly"Beginning Perl for bioinformatics" (which is likely far too basic), which establish some core principles and have a solid reputation in teaching. However, the Perl/Python thing might get confusing (I dunno).

If you are statistics orientated at a guess, O'Reilly's, Bioinformatics data skills... I don't know, but something that would focus on R.

The real goal is core-technical competance because the biological understanding will be taught by your supervisor. Thereafter you might turn your attention to teaming up with a lab for project training.

---
Tags: python
---


# Question
Title: Entrez.efetch returns incomplete genbank records
I am using the biopython Entrez.efetch command to retrieve all features (CDS, mRNA, ...) of some genomes.

In this case (NC\_014649, Acanthamoeba polyphaga mimivirus), it works as expected:

```
from Bio import Entrez, SeqIO
handle = Entrez.efetch(id='NC_014649', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features

[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1181549), strand=1), type='source'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='gene'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2609), strand=1), type='mRNA'),
 SeqFeature(FeatureLocation(ExactPosition(191), ExactPosition(2579), strand=1), type='CDS'),
...

```

However in this case (NC\_010830, Candidatus Amoebophilus asiaticus), it only returns one record for the whole genome.

```
handle = Entrez.efetch(id='NC_010830', db='nuccore', rettype='gb')
records = SeqIO.read(handle, 'gb')
records.features
[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(1884364), strand=1), type='source')]

```

However, if I check the corresponding file in the browser, it does have features: https://www.ncbi.nlm.nih.gov/nuccore/NC\_010830.1

What could be the reason for this inconsistency ? Did I miss something ?

# Answer
> 5 votes
This has been answered before. See: https://stackoverflow.com/a/55402322/6262370

In short, you need to either use `rettype='gbwithparts'` or `rettype='gb', style='withparts'` to download the entire genbank flat file.

# Answer
> 0 votes
I use the `efetch` bash command from **edirect**.

With `efetch -db nuccore -id CP001102.1 -format gb`, I get the *Candidatus Amoebophilus asiaticus 5a2*, complete genome in **GenBank format**.

You might try your command with `CP001102.1` instead of `NC_010830`

---
Tags: genbank, entrez, efetch
---


# Question
Title: What do negative dN or dS mean in cal_dn_ds in Biopython?
I'm noobie in bioinformatics. I use `cal_dn_ds` from `Bio.codonalign.codonseq` with method 'NG86'. Sometimes I get dS=-1. What does it mean? How can number of substitutions be negative? Does it mean that something is wrong with my data?

Upd. 1) Chris\_Rands pointed to sources, I read them already, but the question remain. What is that `ps` variable, and why if it's greater than 3/4, dS = -1?

2) Some more details. `cal_dn_ds` requires sequences from `CodonAlignment` object, which is created by function `Bio.codonalign.build` that requires (quote from docs) `protein alignment and corresponding nucleotides`. I have nucleotide sequences from NCBI at the beginning. So I trim non-coding regions from nucleotide sequences, translate all trimmed sequences, write them to a file, perform multiple alignment of the protein sequences, and then pass the list of aligned protein sequences and the list of nucleotide sequences to the `build` function. I analize protein DKK1, and example sequences that get negatives are *Empidonax traillii XM\_027898752.1* and *Physeter catodon XM\_007110600.1*. But I just want to understand general idea behind dS = -1, what can it mean it different cases.

3) Michael G. gave an interesting suggestion that it means zero, but it's strange because the function returns dN and dS separately, so Biopython might as well return zero, so the user would handle those cases himself. And for some sequences Biopython actually returns normal positive dS and dN = 0 at the same time. If it's an error, what are the possible causes? And why in those error cases Biopython can calculate dN, but dS is -1?

# Answer
dS = -1 .. there are two suggestions

1. I suggest its an error, possibly based on nucleotide calls within the S, check for anomylous nucleotides, particularly over a non-synonymous codon
2. It might mean their are no synonymous mutations. It is not the best idea to assign this a value of 0 in a ratio because it will cause the ratio to fall over when in fact their could be strong selection pressure, ie non-synon mutations. Having an incrementally small value is problematic because dN if positive could anything up to infinity, so -1 might flag this problem.

---

Okay I now know why... You are correct to be concerned with the result because its serious.

My sincere advice to you is do not use any dN/dS ratio where -1 is occurring under any circumstances. On this occasion I'm not going to explain why for two reasons:

1. First and foremost, as a phylogeneticist you should know why, it relates to the fundamentals of the discipline.
2. (Hint) The answer is in my past posts on this forum.

To conclude, I'm pretty nuanced in statistical dogma, you can check my record, IMO most dogma is tribalism because every calculation has weaknesses. HERE HOWEVER, there is a real problem and again I urge you to ditch any data set where -1 is occurring and everything associated with it. That is a fairly dire warning should give you a clue about whats wrong. If you are not prepared to do this, then you must report the result and the situations where it is arising (don't hide it).
> 0 votes

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


# Question
Title: Remove variants that do not map to human genome
*\[This question was also asked on Biostars\]*

I received an hg38 VCF file that's had variants imputed with 1000 genomes. I've encountered some issues with the VCF; REF alleles that do not align to a reference genome, ALT alleles that do not appear to be reported anywhere in the literature, and, most recently, variants that flat-out do not align to the human genome (variants on chr19 with bp-pos 100 million+ when the whole chromosome is in the 50 million bp range).

I've worked out hack-y solutions to most of the issues that I've encountered, but this latest one has been an issue for me. I only detected these variants when I ran VEP and it flagged them as not mapping to the genome. As such, I'm more or less removing these variants one at a time using `grep -v`. I'd like a solution where I can just remove any variants from the vcf that appear to map to regions that do not exist in the human genome. Bonus points if the solution also encompasses some of the other issues I mentioned, although I think I've already found solutions to those. Is there anything out there that does this?

# Answer
Answer I received on Biostars that appears to work as expected:

`bcftools norm -Ob --check-ref wx -f <reference_fasta> <vcf.gz> -o <output_name>`
> 1 votes

---
Tags: vcf, variant-calling, variation, quality-control
---


# Question
Title: bowtie2 options when mapping stranded single end reads
I'm trying to map RNA-Seq reads generated using the `NEB Ultra Directional kit is first strand reversed` protocal. I'm using bowtie-2.3.5\_1 to align against the mus musculuous (mm10) genome. The reads are single ended.

It seems that `--fr` and `--rf` are relevant only for paired end options. I tried these options as well as not specifying either (default).

E.g.

```
$ bowtie2 --fr \
        --threads 20 \
        --very-sensitive \
        -x /reference/mus_musculus/mm10/ucsc_assmebly/illumina_download/Sequence/Bowtie2Index/genome \
        -U trimmed_out/my_trimmed.fq.gz \
        -S bowtie_fr_1.sam

$ samtools sort -n -o bowtie_fr_sort_1.sam bowtie_def_1.sam

```

I then compared the output of `bowtie_rf_sort.sam` `bowtie_default_sort.sam` and `bowtie_fr_sort.sam`. They are basically identical excluding the header.

**Question**:

Does strandedness matter in this alignment? Perhaps I should care about the `--nofw/--norc` options? Am I missing something here?

# Answer
> 1 votes
I repeat what user172818 says...don't waste your time on bowtie or tophat. Use STAR or HISAT2, or a pseudoaligner like kallisto with a transcriptome.

---
Tags: rna-seq, read-mapping, bowtie2
---


# Question
Title: Why are my kallisto and salmon results differing so much just for lncRNA transcripts?
I am running some analysis on an RNA-seq dataset. I have a list of transcripts that are potential lncRNA for which I ran both Kallisto and Salmon aligners. The input data for index building and quantification includes mRNA as well as these potential lncRNA.

Commands for one replicate:

Salmon

```
salmon index --keepDuplicates -t transcripts.fasta -i index_directory/ -k 31
salmon quant --seqBias --gcBias -i index_directory/ --libType ISR --auxDir index_directory/ -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz -o output_name

```

Kallisto

```
kallisto index -i kallisto_index.idx transcripts.fasta
kallisto quant -i kallisto_index.idx -o output_name --rf-stranded sample_R1.fastq.gz sample_R2.fastq.gz

```

Then I used the tpm values, averaged replicates (3) and plotted the number of transcripts per condition that had a tpm \> 2. Same result for tpm \> 5.

The following graph shows the for mRNA:

And for lncRNA:

Again, the analysis were done with both lncRNA and mRNA transcripts together in the same fasta file.

The overlap is not shown here but is always smaller than the count for Kallisto, for the lncRNA.

In literature both tools seem to have similar results. I cannot find a reason for this. If there would be an error somewhere I would expect to see also differences in the conding RNAs which does not happen.

# Answer
There has been a detailed comparison of kallisto and salmon for lncRNA transcripts by Zheng et al https://www.biorxiv.org/content/biorxiv/early/2018/01/09/241869.full.pdf

It shows them to be near identical. This suggests there is some problem with the way you made the comparison/plot.
> 2 votes

---
Tags: rna-seq, salmon, kallisto
---


title: blast command line

text:
# Question
Title: blast command line
I would to use blastx from command line to analize a library produced with RepeatScout in format .fasta to find new transposable elements. I use this command `mylaptop$ tblastx OB150DEC_2_repeats.fa`

but it reports an error:

```
USAGE
  tblastx [-h] [-help] [-import_search_strategy filename]
    [-export_search_strategy filename] [-db database_name]
    [-dbsize num_letters] [-gilist filename] [-seqidlist filename]
    [-negative_gilist filename] [-entrez_query entrez_query]
    [-db_soft_mask filtering_algorithm] [-db_hard_mask filtering_algorithm]
    [-subject subject_input_file] [-subject_loc range] [-query input_file]
    [-out output_file] [-evalue evalue] [-word_size int_value]
    [-qcov_hsp_perc float_value] [-max_hsps int_value]
    [-xdrop_ungap float_value] [-searchsp int_value] [-sum_stats bool_value]
    [-max_intron_length length] [-seg SEG_options]
    [-soft_masking soft_masking] [-matrix matrix_name]
    [-threshold float_value] [-culling_limit int_value]
    [-best_hit_overhang float_value] [-best_hit_score_edge float_value]
    [-window_size int_value] [-lcase_masking] [-query_loc range]
    [-strand strand] [-parse_deflines] [-query_gencode int_value]
    [-db_gencode int_value] [-outfmt format] [-show_gis]
    [-num_descriptions int_value] [-num_alignments int_value]
    [-line_length line_length] [-html] [-max_target_seqs num_sequences]
    [-num_threads int_value] [-remote] [-version]

DESCRIPTION
   Translated Query-Translated Subject BLAST 2.5.0+

Use '-help' to print detailed descriptions of command line arguments
========================================================================

Error: Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: OB150DEC_2_repeats.fa

```

I don't understand what does it mean `too many positional arguments`

# Answer
> 1 votes
A positional argument is one that doesn't follow something like `-foo` or `--bar`. That is, the variable that it should internally map to depends on it's position in the command. You're looking for something like `tblastx -db something -query OB150DEC_2_repeats.fa`

---
Tags: blast
---


# Question
Title: Align legends to each heatmap in a ComplexHeatmap with 2 heatmaps
I have created a ComplexHeatmap containing 2 `Heatmap`s, each with their own legend. The legend on the first Heatmap is continuous, whereas that of the second is discrete.

When I concatenate these `Heatmap`s vertically, and `draw` the `HeatmapList`, their legends are auto-aligned to the center of the plot. Is there any way I can have these legends positioned centered to the respective heatmaps and not to the center of the plot?

This is what I have:

Here's what I'd like to have:

How can I go about solving this problem? I'd really like to avoid creating customized `Legend`s if I can, but if that's what it takes, that works too.

Thank you in advance for any help/pointers!

Note: I tried adding the `ComplexHeatmap` tag, but it doesn't exist yet and I don't have the reputation to create new tags yet.

Edit: Since this is a Bioconductor package, I also created an identical post on the Bioconductor support forum

Dummy data and code:

```
devtools::install_github('jokergoo/ComplexHeatmap');
library(ComplexHeatmap);
library(viridis);

dummy_top_mat <- structure(c(0.4563, 0.2211, 1.2235, 0.067, 1.6859, 1.0936, 0.4533,
1.6844, 1.0039, 0.5402, 0.6841, 1.498, 1.042, 0.1711, 0.3565,
2.2814, 3.0516, 1.5012, 0.5367, 3.0846, 1.1909, 0.235, 1.4266,
0.4858, 2.9054, 0.3733, 0.6902, 0.8555, 0.4234, 2.6778, 0.1568,
0.5556, 2.0172, 0.8034, 2.2897, 0.1166, 3.8033, 0.1431, 2.0606,
1.2725, 1.5365, 0.4123, 1.2087, 1.1264, 0.8334, 1.1943, 1.58,
1.5849, 0.3004, 0.3722, 0.0362, 0.0532, 1.4867, 0.4053, 0.3615,
0.0897, 1.3217, 1.1447, 1.3058, 0.1903, 0.1067, 0.9482, 1.3382,
3.2955, 0.391, 1.0418, 0.2041, 1.208, 1.5857, 3.5313, 0.472,
1.389, 0.2143, 0.0226, 0.029, 0.444, 2.0521, 0.3955, 0.3495,
0.5062, 1.3236, 1.3234, 0.7111, 0.1176, 2.2223, 1.2073, 0.3964,
2.1175, 0.3382, 0.2816, 0.71, 3.1417, 0.2402, 0.5793, 0.7662,
1.6782, 0.0986, 0.087, 0.5447, 2.6672, 1.2498, 1.0676, 1.8608,
1.8146, 0.1422, 0.4221, 0.0303, 0.9541, 0.7358, 1.7664, 1.5144,
0.2034, 0.9366, 0.7837, 0.3284, 0.1477, 1.8306, 1.3564, 0.1126,
0.0171, 2.9858, 0.0233, 0.2796, 0.6995, 1.6081, 0.215, 1.7093,
0.5178, 1.7061, 2.473, 1.8912, 0.7661, 4.4102, 1.2963, 0.6542,
0.4281, 0.4491, 0.6, 0.4076, 1.6869, 0.4747, 3.9823, 1.1226,
2.7355, 2.7036, 0.2241, 0.983, 1.0992, 1.4736, 0.1584, 0.2995,
0.2272, 0.5744, 0.9314, 0.6924, 0.0812, 1.361, 0.727, 0.1525,
1.3367, 0.566, 2.7801, 0.2349, 3.2655, 1.0675, 0.449, 0.5411,
0.8291, 0.52, 0.0507, 0.6538, 0.4636, 1.2063, 0.9784, 0.1925,
0.0756, 0.136, 1.2529, 0.317, 0.0281, 0.8668, 3.4138, 5.3898,
0.521, 0.087, 3.4819, 0.1114, 0.0061, 0.9804, 1.139, 1.4159,
1.0297, 0.6503, 1.0828, 2.3527, 0.057, 0.5602, 0.4017, 0.9985,
0.8673), .Dim = c(10L, 20L), .Dimnames = list(c("gene01", "gene02",
"gene03", "gene04", "gene05", "gene06", "gene07", "gene08", "gene09",
"gene10"), c("sample01", "sample02", "sample03", "sample04",
"sample05", "sample06", "sample07", "sample08", "sample09", "sample10",
"sample11", "sample12", "sample13", "sample14", "sample15", "sample16",
"sample17", "sample18", "sample19", "sample20")));

dummy_bottom_mat <- structure(c("Positive", "Positive", "Positive", "Positive", "Negative",
"Positive", "Positive", "Positive", "Positive", "Negative", "Positive",
"Positive", "Positive", "Positive", "Positive", "Negative", "Positive",
"Negative", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Negative", "Positive", "Positive", "Positive", "Positive", "Positive",
"Negative", "Positive", "Positive", "Positive", "Negative", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Negative",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Negative", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive", "Positive", "Positive", "Positive", "Positive", "Positive",
"Positive"), .Dim = c(6L, 20L), .Dimnames = list(c("marker1",
"marker2", "marker3", "marker4", "marker5", "marker6"), c("sample01",
"sample02", "sample03", "sample04", "sample05", "sample06", "sample07",
"sample08", "sample09", "sample10", "sample11", "sample12", "sample13",
"sample14", "sample15", "sample16", "sample17", "sample18", "sample19",
"sample20")));

hmap1 <- Heatmap(dummy_top_mat, col = colorRampPalette(viridis(15))(100), name = 'log2 (TPM)', na_col = 'grey60', cluster_rows = FALSE, cluster_columns = FALSE, row_names_side = 'left', column_title_side = 'top', rect_gp = gpar(col='white', lwd=0.5), height=unit(10*0.75,'cm'), width = unit(20*0.75, 'cm'), column_title = 'Gene Expression');

hmap2 <- Heatmap(dummy_bottom_mat, col = c('Negative'='red', 'Positive'='blue'), na_col = 'grey60', row_split = factor( c('gp1','gp1','gp1','gp1','gp2','gp2'), c('gp1','gp2'), ordered = TRUE), name = 'Marker\nStatus', rect_gp = gpar(col='white',lwd=1), row_title = NULL, height = unit(6*0.75,'cm'), width = unit(20*0.75, 'cm'), row_names_side = 'left', column_title = 'Biomarkers');

draw(hmap1 %v% hmap2, ht_gap = unit(1, 'cm'), auto_adjust = FALSE);

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS  10.14.4

Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Users/ram/miniconda3/lib/R/lib/libRblas.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] viridis_0.5.1        viridisLite_0.3.0    ComplexHeatmap_2.1.0 forcats_0.4.0
 [5] stringr_1.4.0        dplyr_0.8.0.1        purrr_0.2.5          readr_1.3.1
 [9] tidyr_0.8.3          tibble_2.1.1         ggplot2_3.1.1        tidyverse_1.2.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.0          lubridate_1.7.4     lattice_0.20-38     circlize_0.4.6      prettyunits_1.0.2
 [6] png_0.1-7           ps_1.3.0            rprojroot_1.3-2     assertthat_0.2.1    digest_0.6.18
[11] R6_2.3.0            cellranger_1.1.0    plyr_1.8.4          backports_1.1.4     reprex_0.2.1
[16] httr_1.4.0          pillar_1.3.1        GlobalOptions_0.1.0 rlang_0.3.1         lazyeval_0.2.2
[21] readxl_1.3.1        rstudioapi_0.10     callr_3.2.0         GetoptLong_0.1.7    desc_1.2.0
[26] devtools_2.0.2      munsell_0.5.0       broom_0.5.2         compiler_3.5.1      modelr_0.1.4
[31] pkgconfig_2.0.2     pkgbuild_1.0.3      shape_1.4.4         tidyselect_0.2.5    gridExtra_2.3
[36] crayon_1.3.4        withr_2.1.2         nlme_3.1-137        jsonlite_1.6        gtable_0.3.0
[41] magrittr_1.5        scales_1.0.0        cli_1.1.0           stringi_1.4.3       remotes_2.0.4
[46] fs_1.3.1            xml2_1.2.0          generics_0.0.2      rjson_0.2.20        RColorBrewer_1.1-2
[51] tools_3.5.1         glue_1.3.0          hms_0.4.2           processx_3.3.1      pkgload_1.0.2
[56] parallel_3.5.1      clue_0.3-57         colorspace_1.4-1    cluster_2.0.7-1     sessioninfo_1.1.1
[61] rvest_0.3.3         memoise_1.1.0       haven_2.1.0         usethis_1.5.0

```

# Answer
> 1 votes
Turns out, this is not possible. Concatenating two heatmaps causes their legends to be grouped, and while they can be positioned together at the top, bottom, at the left or the right sides, they cannot be separated spatially.

The discussion can be viewed at the bioconductor support link from the question. Relevant snippets:

> No, since you put the two heatmaps as one plot, the legends are grouped. You can either edit in Inkscape or make two heatmaps in two plots. You can get the column (order) of the first heatmap ... and assign to the second heatmap

> The thing is if you add two heatmaps, the legends are always plotted together. One way I can think to improve this is to move the legends to the bottom of the heatmap.

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


title: Multi-line efetch digging

text:
# Question
Title: Multi-line efetch digging
I'm using efetch to dig out a middle "envelope" protein within a viral genome specifically within fasta format. The code below works but may not be robust for an entire database. Thus I want to replace the two "perl" lines with formal code within the Entrez package. I tried using `efilter -query` but couldn't get it to work. Genome positions are one approach but are not absolute from one virus to another. Rewriting this may be the best option.

```
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" |\
xtract -pattern INSDSeq -element INSDQualifier_value |\
perl -p -e 's/^(.+?)\s+taxon.*.*\s+\b([A-Z]{300}[A-Z]+).*/>$1\n$2/'|\
perl -p -e "s/[\t \;\/]+/_/g"

```

Output

```
>InfluenzaAvirus_A_PuertoRico_8_1934
MDVNPTLLFLKVPAQNAISTTFPYTGD.....

```

I will provide a more detailed explanation shortly. Any ideas welcome.

# Answer
Not exactly what you were looking for but may be you can use something like this?

```
ncbi2=NC_002021
efetch -db protein -format gpc -id "$ncbi2" \
    | xtract -pattern INSDSeq_feature-table \
        -group INSDFeature \
        -if INSDFeature_key -equals 'mat_peptide' \
        -tab '\n' \
        -element INSDFeature_location,INSDQualifier_value \
    | grep 'Segment 2'

```

For explorations of the INSDSet XMLs, `xtract` has `-insd` option that you can use. In this particular case, you can do something like this:

```
efetch -db nuccore -id "$ncbi2" -format gpc \
    | xtract -insd mat_peptide product peptide

```

For more information, checkout EDirect documentation and some example queries
> 2 votes

---
Tags: covid-19, bash, sars-cov-2, perl, efetch
---


title: Tabix queries spanning chromosomes

text:
# Question
Title: Tabix queries spanning chromosomes
I can do something like this:

```
$ tabix ${url} ${range}

```

where `${range}` spans a region on a single chromosome, *e.g.*, `chrA:X-Y`.

Can a `tabix` query span multiple chromosomes, *e.g.*, where `${range}` looks like `chrA:X-chrN:Y` (or similar), which returns data between `chrA:X` and `chrN:Y` (using either a lexicographical or natural sort order between `chrA` and `chrN`)?

*Note: I do not want to pass in a BED file of intervals. If it can't be done without passing in a separate file, that's fine. I only want to provide a one-liner to give to someone who can run `tabix` on their end, to do this type of query directly and remotely.*

# Answer
> 3 votes
Looking at the htslib source code, it seems like the answer is no, a single region specifier cannot span multliple chromosomes.

Seems like the closest you can get is to enumerate all the chromosomes as separate regions on the command line. So to get everything between `chrA:X` and `chrN:Y`, you would have to write:

```
tabix inputfile.gz chrA:X chrB chrC ... chrN:1-Y

```

(If you omit the end coordinate like `chrA:X`, the region encompasses `chrA` from the specified coordinate to the end. If only specify the chromosome name without the coordinates, like `chrB`, the region encompasses the whole named chromosome, and for the last chromosome, you can specify up to the desired point by starting from the beginning: `chrN:1-Y`.)

---
Tags: tabix
---


# Question
Title: Which values from Seurat::FetchData function are to be used for correlation analysis between genes?
I want to perform a correlation test between genes in on my single cell RNA seq data set. I perfomed the differential expression analysis using the **Seurat** version 2 package, after performing stages of normalisation, scaling, PCA, TSNE analyses and clustering. My correlation plot is far from what I expect and I want to know if the expression values obtained with **Seurats** V2's **FetchData** function (imputed, scaled, or raw) may have an effect on correlation analysis' outcome. This is what I did:

```
# load Seurat v2
library(Seurat)

# load corrplot package
library(corrplot)

# obtain expression values of the genes I want to test rom the Seurat object
expression_matrix <- FetchData(object, vars.all = c("gene1", "gene2", "gene3"))

expression_matrix <- round(expression_matrix , 1)

cor_expression_matrix <- cor(expression_matrix)

corrplot(cor_expression_matrix , is.corr = T, method = "square")

```

Which gene expression values should be best used for the analysis of correlation between genes: imputed, scaled, or raw?

Thank you in advance for your kind response.

# Answer
> 0 votes
You should use normalized data. The function has a parameter slot="data" that let's you select to access the normalized data.

Any reason for rounding to 1 decimal place here? I imagine that would have a substantial impact on the correlation you compute.

---
Tags: seurat, single-cell, correlation
---


# Question
Title: Alternative splicing positions on protein
I have mass-spectrometry data (Human cellsanalyzed and annotated by MaxQuant). I want to find which peptides are covering positions of former splice junctions (on the mRNA sequence), based on protein id and peptide positions.

Where could I find an appropriate data? I want to know, based on peptide-sequence, annotaion and position on protein, whether this peptide is a single exon, or if it contained areas of former splice junctions.

Itai

# Answer
If you are considering the human genome as "data" (not mRNA), I would run your single peptide via tblastn which will take your protein(s) translate them to nucleotides and run a blast against Genbank. The splice sites will be immediately apparent because they force indels in your query. If there are no indels there are no introns (splice sites). If you want to process several proteins, one solution is to use eblast, a command line version of tblastn (Anaconda package). You can confirm this against a full genome annotation (note, I don't work on humans).

Alot of work has been done on splice sites, a friend worked on this, but I think you just want to focus on a restricted number of target peptides in your mass spec. output.

mRNA should be avoided IMO if you are looking for splice sites, i.e. intron boundaries.
> 0 votes

---
Tags: proteins, isoform
---


# Question
Title: Differential expression between two types of cells of specific marker gene expression
I want to define the identity of some cells based on a gene marker (for example the colour yellow\< EYFP gene). One group of cells has the EYFP gene the other doesn't and I want to find the differential gene expression between those groups that are found on the same cell sample. Do you have any tips for a code in Seurat? Thanks

# Answer
> 1 votes
Take a look at the section **"Finding Differentialy Expressed Feautres (cluster biomarkers)"** of Seurat tutorial.

You have to define your clusters (if not computed automatically) by modifying the ***ident*** slot of the Seurat object.

Then you call the function *FindMarkers* of your Seurat object (to define markers of a given cluster against the rest) or *FindAllMarkers* (to look for all the combinations)

---
Tags: seurat
---


# Question
Title: How to merge .fastq.qz files into a single .fastq.gz with their same id without losing any content in parallel
I have a large number of .fastq.gz files of different lane and reads. I have to merge them each reads group files into single .fastq.gz files.

\**eg:

**1st type**
NA24694\_GCCAAT\_L001\_R1\_001.fastq.gz
NA24694\_GCCAAT\_L001\_R1\_002.fastq.gz
…....
…....
NA24694\_GCCAAT\_L001\_R1\_033.fastq.gz

**2nd type**
NA24694\_GCCAAT\_L001\_R2\_001.fastq.gz
NA24694\_GCCAAT\_L001\_R2\_002.fastq.gz
…....
…....
NA24694\_GCCAAT\_L001\_R2\_027.fastq.gz

**3rd type**
NA24694\_GCCAAT\_L002\_R1\_001.fastq.g
NA24694\_GCCAAT\_L002\_R1\_002.fastq.g
…....
…....
NA24694\_GCCAAT\_L002\_R1\_040.fastq.g

**4th type**
NA24694\_GCCAAT\_L002\_R2\_001.fastq.gz
NA24694\_GCCAAT\_L002\_R2\_002.fastq.gz
…....
…....
NA24694\_GCCAAT\_L002\_R2\_040.fastq.gz

so now i need to merge all the corresponding files into single file.

**Output:**
EA00694\_GCCAAT\_L001\_R1.fastq.gz
EA00694\_GCCAAT\_L001\_R2.fastq.gz
EA00694\_GCCAAT\_L002\_R1.fastq.gz
EA00694\_GCCAAT\_L002\_R2.fastq.gz

I tried ""cat'" but its not possible to do for large files. Could any one help me to run in script parallel in perl or python or shell. or any linux command which can be used for large file.

Thanks all.

# Answer
> 10 votes
All you need is `cat`. You won't find any better tool for a simple job like this. Just run:

```
cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz
cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz
cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz
cat NA24694_GCCAAT_L001_R4*fastq.gz > EA00694_GCCAAT_L002_R2.fastq.gz

```

That should work just fine. You say you tried `cat` and it "didn't work" but since you don't tell us how it failed, I can't really help. The only issues I can think of is that either there are too many files (as in several hundred thousand, whatever the value returned by `getconf ARG_MAX` on your system which is 2097152 on mine) or, more likely, you are running out of disk space.

If it's a disk space issue, you might be able to get around it by adding each file and then deleting it:

```
for i in 1 2 3 4; do
    for file in NA24694_GCCAAT_L001_R${i}_*fastq.gz; do
        cat "$file" >> EA00694_GCCAAT_L001_R${i}.fastq.gz && rm "$file"
    done
done

```

***IMPORTANT***: the command above is destructive. It will delete each file after it has been added to the new one. If everything works fine that's not a problem, but if I've made a mistake or your file names are slightly different you might lose your data. So I strongly urge you to simply run the `cat` commands in the first section on a machine with enough disk space to deal with them.

# Answer
> 3 votes
We sequence DNA on daily basis. Like `terdon` said, you can use `cat` to group all the FastQ of all your sample whether they are bgzipped or not like so:

```
root@slurm1: ~  # ls
am224_S6_L001_R1_001.fastq.gz
am224_S6_L001_R2_002.fastq.gz
am224_S6_L002_R1_003.fastq.gz
am224_S6_L002_R2_004.fastq.gz
am232_S10_L001_R1_001.fastq.gz
am232_S10_L001_R2_002.fastq.gz
am232_S10_L002_R1_003.fastq.gz
am232_S10_L002_R2_004.fastq.gz

```

You could do:

```
# initiate am empty array
listSamples=()

# loop to add the unique samples to an array
for i in *.fastq.gz;
do
    if [[ "$sampleName" =~ "${listSamples[@]}" ]]; then
        sampleName=$(echo $i | cut -d'_' -f1)
        listSamples+=("$sampleName")
    fi
done

# loop in the newly created array of samples and group them by reads number
for i in "${listSamples[@]}"
do
    cat "$i"_*_R1_* > "$i"_1.fastq.gz
    cat "$i"_*_R2_* > "$i"_2.fastq.gz
done

```

However, I highly suggest you to compress your files as `.bz2` to save storage as much as possible. You can achieve this using the multithreaded tool lbzip2. Assuming that you have 10 CPUs available, you can do:

```
for i in "${listSamples[@]}"
do
    bzcat "$i"_*_R1_* | lbzip2 -n 10 --best > "$i"_1.fastq.bz2
    bzcat "$i"_*_R2_* | lbzip2 -n 10 --best > "$i"_2.fastq.bz2
done

```

Hope this helps !

---
Tags: bash, shell, data-management
---


# Question
Title: Proper use of BWA MEM on multiplexed GBS sample
I have a multiplexed lane of GBS sequencing reads as a fastq file. I understand the first step is to demultiplex and trim the adapter sequences from the reads.

This yields many individual fastq files that correspond to individual libraries.

I would like to genotype my samples, but I need to map the reads first. My question is how should I perform the BWA MEM alignment?

Do I perform BWA MEM on each of the resulting fastq files that correspond to individual libraries...

OR

I should I perform BWA MEM on the large original fastq file that contains all my sequence reads before demux and read trimming?

Forgive my ignorance, I've only ever used TASSEL (which works great) and I'm trying to expand my understanding and flexibility.

# Answer
> 1 votes
The first thing to do is to group the FastQ files together following this post. Once you have your reads grouped, you could trim your files before mapping. Then you could perform your BWA alignment. Since you mentioned `MEM`, I suppose your read length should be between 70-100bp as suggested by the author:

```
It consists of three algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
The first algorithm is designed for Illumina sequence reads up to 100bp, while the rest two for
longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar features such
as long-read support and split alignment, but BWA-MEM, which is the latest,
is generally recommended for high-quality queries as it is faster and more accurate.
BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.

```

Then you could run your BWA MEM, supposing you have 10 CPUs available, like so:

```
bwa mem \                                  # bwa mem algorithm
-M \                                       # mark shorter split hits as secondary (for Picard compatibility).
-t 10 \                                    # CPUs thread
reference.fasta \                          # reference genome to be aligned to
< sample_r1.fastq \                        # fastQ read 1 (can be .gz)
< sample_r2.fastq | \                      # optional FastQ read 2 (can be .gz)
samtools view -bS -@ 10 > sample.bam       # pipe to samtools to convert SAM to BAM on the fly.

```

Hope this helps !

# Answer
> 0 votes
Trimming is something you always want to do before mapping to make sure that the bases at each position are of sufficient quality, and to make sure that you are not aligning your adapters to the reference genome.

I am not entirely sure what the biological difference is represented in your individual libraries (that depends on your research question). But most likely you want to demultiplex the reads before mapping with BWA as well. As a result you would indeed get a sam/bam file per library. Let me know if anything is still unclear!

---
Tags: bwa, genotyping
---