# Question
Title: Why do I get so many insertions from Minimap2 on my Nanopore WGS?
I'm a starting my analysis on nanopore whole genome sequencing. I start my analysis from this popular Github.

The sample I downloaded was WGS for NA12878, so I would assume it's alignment to GRCh38 shouldn't be that bad.

But ... I'm getting lot's of insertions...

I'm not sure it was caused by the sequencing platform, the MiniMap aligner or my bugs.

I supply all my commands, starting from sample download.

```
wget https://s3.amazonaws.com/nanopore-human-wgs/rel5-guppy-0.3.0-chunk10k.fastq.gz
gunzip rel5-guppy-0.3.0-chunk10k.fastq.gz
wget http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -d GRCh38_full_analysis_set_plus_decoy_hla.fa.index GRCh38_full_analysis_set_plus_decoy_hla.fa
minimap2 -ax map-ont -t 20 GRCh38_full_analysis_set_plus_decoy_hla.fa.index rel5-guppy-0.3.0-chunk10k.fastq | samtools view -bS - | samtools sort -@ 5 > rel5-guppy-0.3.0-chunk10k.bam
samtools index rel5-guppy-0.3.0-chunk10k.bam

```

# Answer
Insertions and deletions are the dominant error mode of long read sequencing, including nanopore sequencing. What you see is not unexpected. Things may have improved by now if you would download the raw fast5 data and repeat the basecalling. There is no need to `gunzip` the fastq.gz prior to alignment. Your commands for alignment look alright, except that (if you are using a recent version of samtools) you don't really need samtools view and can pipe directly to samtools sort, using -o to specify the output format:

```
minimap2 -ax map-ont -t 20 <index> <fastq> | samtools sort -@5 -o alignment.bam

```

In a recent version of igv you can right-click the alignment and chose a "Quick consensus mode" (or something like that), which will make your alignments at least a bit prettier. Note how the consensus sequence (top of the graph) is quite clean - while individual reads may be noisy this is mostly randomly distributed and as such ironed out in consensus.
> 8 votes

---
Tags: sequence-alignment, ngs, nanopore, errors
---


# Question
Title: RNA-seq gene/transcript read counts database for Mouse
Is there any RNAseq gene/transcript read count database for the mouse? I already know about ARCHS4, looking for some other source.

Thank you

# Answer
EBI has an expression atlas, including mouse: https://www.ebi.ac.uk/gxa/experiments?species=mus+musculus

If you click on one of the individual experiments you can download tables of TPM or RPKM per gene.

There are currently about 1000 individual experiments in different tissues etc.
> 2 votes

---
Tags: rna-seq, reads, data-mining
---


# Question
Title: Recommendation for a binary dataset for Computational Biology open access (as in a dataset with 1 and 0 that can be used to apply ML techniques?
I have a couple of modified ML methods that I want to test with biological binary data - specifically binary gene expression data - e.g. https://eprints.soton.ac.uk/69359/1/Thesis.pdf that might include a variable that could be used to be the target variable in ML classification. What I mean by binary data is the features are all either 0 or 1, so I expect a large dataset e.g this thesis. I am looking for any suggestions.

# Answer
> 2 votes
Kaggle has a cats v dogs dataset if you want to try your hand at biological image classification: https://www.kaggle.com/c/dogs-vs-cats

The DREAM challenge competitions are usually specifically oriented to solving biological classification or learning problems with ML: http://dreamchallenges.org/

---
Tags: machine-learning, public-dataset
---


# Question
Title: Selecting 65000 SNPs where AF is close to 0.5 in all or most populations
I am evaluating the tool somalier (https://github.com/brentp/somalier) and I need to create a list of about 65,000 SNPs where the allele frequency (AF) is as close to 0.5 as possible across the most representative set of populations possible with today's data.

My alignments are in GRCh38, so I need a vcf input file for somalier also in GRCh38, not the one that is available in the somalier website, which is in GRCh37.

I am starting as input with the vcfs from the 1000genomes project in GRCh38 positions, and I would like to come up with a filtering strategy that: (1) Sorts all SNPs by the AF being close to 0.5 in as many populations as possible. (2) Takes the first 65,000 SNPs from that list and produces a vcf file from them.

My simplistic strategy so far has been to take the global AF value if it's between 0.3-0.7 and print it if it's a SNP with this perl one-liner:

```
gunzip -c my.vcf.gz | grep -v '#' | perl -lne 'if ($_ =~ /\;AF=0\.[34567]/) { if ($_ =~ /\w+\t\d+\t\S+\t\w\t\w\t/) {print $_} }'

```

Any recommendations?

# Answer
> 5 votes
Just use bcftools view for filtering:

```
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf

```

To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use `head -n`.

```
$ bcftools view -h input.vcf.gz|wc -l
255

$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | head -n 65255 > output.vcf

```

# Answer
> 3 votes
Using `bcftools` is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single `zgrep` call:

```
zgrep -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

```

And to limit to the first 65000 results:

```
zgrep -m 65000 -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz

```

Finally, to keep the header of the original vcf unchanged, you could do:

```
zgrep -m 65000 -P '^(#|[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7])' my.vcf.gz

```

I also suggest you use hg19 instead of hg38 data since the best available source for frequencies, the gnomAD database, is based on hg19 and the tool you want to test also provides a site file for hg19. But that may not be possible for you.

If your VCF file has multiple AF entries for different sub-populations and you want to use those, you'll probably need to script something. We can't really help with that unless you show us some example input lines and start writing the basic script yourself though.

---
Tags: vcf, snp, genomics, variation
---


title: Calculate average parental genotype

text:
# Question
Title: Calculate average parental genotype
I am interested in creating an averaged value for genotype across two parents:

E.g. Input (three unrelated couples, 1, 2 and 3):

```
Individual    SNP1    SNP2    SNP3    SNP4 ...
Father_1      1       0       1       2    ...
Mother_1      0       0       2       1    ...
Father_2      1       0       1       1    ...
Mother_2      2       1       0       1    ...
Father_3      1       2       0       1    ...
Mother_3      1       .       0       2    ...

```

E.g. Output (the midparental genotypes for couple 1, 2 and 3):

```
Individual    SNP1    SNP2    SNP3    SNP4 ...
Midparent_1   0.5     0       1.5     1.5  ...
Midparent_2   1.5     0.5     0.5     1    ...
Midparent_3   1       .       0       1.5  ...

```

I can think of a way to do this using e.g. a PLINK .raw file and manipulating that in R, but am concerned that would become memory- and storage-heavy when assessing 100000s of SNPs and 1000s of pairs.

Does anyone know of an existing tool that can do this?

Thanks

# Answer
Beside the correct question by @Kamil, this sounds like a job for datamash.

First let's get the number of columns in the file

```
$ head -n1 genotypes.txt| awk -v FS="\t" '{print NF}'

```

Now lets calculate the average per column with `datamash`

```
$ datamash -H mean 2-<n> < genotypes.txt > mean.txt

```

Replace `<n>` with the number of columns we received with command above.

For your example the output will look like this:

```
mean(SNP1)  mean(SNP2)  mean(SNP3)  mean(SNP4)
0,5 0   1,5 1,5

```

One can remove the `mean` with `sed`

```
$ sed -i -e 's/mean(//g' -e 's/)//g' mean.txt

```
> 2 votes

# Answer
These are unusual parents because they all appear haploid. Anyway, if I was doing this I would calculate the allele frequency per progeny such as makefreq in adegenet in R-cran (just my personal thoughts).

**Problem** The issue is that you are using a continuous distribution 1,2,3 etc... for catagorical data and will result in confusing homozygotes with heterozygotes.

***For example***

* If you have allele 3 (father) and allele 1 (mother) the "progeny" wil be "allele" 2
* If you have allele 2 (father) and allele 2 (mother) the "progeny" will be "allele" 2

However, one is homozygous and the other is heterozygous and the allelic background for each is quite differnt.

The alternative is to use letters rather than numbers. You can always use "OneHotKeyEncoder" to convert back from AA, aa etc... to numerical catagories for analysis later (its a standard feature in Python).

If you used letters rather and then converted to numbers every "progeny" would have its own numerical catagory (but you would lose homozygous vs. heterozygous classification). It is standard practice to switch from letters to numbers and the "encoder" will ensure no mistakes are made. The other way around this is to designate the variable as a string, rather than as an integer... this might be more to your liking, but you are stuck in base 10.

With regards the throughput, 100000s 1000 individuals, a 100 million SNPs is doable via R, the operations are simple.

Just to clarify: it would be easy to script to a solution (even from base code), but I'm not sure about how to interpret the output. All I'm suggesting is you have another think about the goal of this analysis.
> 1 votes

---
Tags: snp, plink, data-management
---


# Question
Title: RSeQC provided bedfile does not work with read_distribution.py
I frequently make use of RSeQC, I used to use it with my own bed file, generated using gtf2bed, but now I thought I'd use the one provided by RSeQC themselves from here, I got the file hg38\_Gencode\_V28.bed.gz, when looking at it with IGV it looks similar to my self generated bedfile however, read\_distribution.py (a script in RSeQC), does not report any reads in any of the categories (CDS\_exons, Introns, etc).

Any idea why?

Edit: I see that my own bed file (from the GTF) had ENST identifiers like this ENST00000466557 whereas the file from RSeQC supplies them like ENST00000611746.1 (with decimals) that seems to be the only difference as far as I can see.

# Answer
> 2 votes
Ok, I figured it out already, it's and old issue I come accross every now and then. The original bed file lists the chromosomes as chr1,chr2,...,chrX whereas my bam files uses 1,2,...,X

I fixed the bed file like this (python3):

```
import pandas as pd
bedFile = '/hg38_Gencode_V28.bed'
bedFile = pd.read_csv(bedFile, sep='\t', header=None)
bedFile[0] = [value.strip('chr') for value in bedFile[0]]
newBedFile = 'hg38_Gencode_V28_reformatted.bed'
bedFile.to_csv(newBedFile, sep='\t', header=None, index=None)

```

And now read\_distribution.py works.

---
Tags: rna-seq, qc, rseqc
---


# Question
Title: How to model the effect of acid/base on pH
I would like to model a biological compartment with sysBio R package. It is fine but I'm having a difficult time with modeling the change in pH due to acidosis. I tried several things like:

```
addMAreaction(glu, "pH = null", "acidosis", "base excess")
addMAreactRate(glu, "acidosis", "assigned", "ka*pH")
addMAreactRate(glu, "base excess", "assigned", "ka*pH*-1")

```

which caused continuously decreasing pH levels. I'll leave out the intermediate steps but the simplified version of the last code looks like this:

```
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "pH", 7.4)
addSpecies(glu, "Acid", 0)

addMAreaction(glu, "Glucose = Lactate + 0.01*Acid", "forwardL",
              "backwardL", name = "glycolysis, sort of")

addMAreactRate(glu, "forwardL", "assigned", "kgl*Glucose")
addMAreactRate(glu, "backwardL", "assigned", "ka*pH+klg*Lactate")

addParameters(glu, "kgl", 0.01)
addParameters(glu, "klg", 0.01)
addParameters(glu, "ka", 0.05)
addParameters(glu, "a", 0.001)

addRule(glu, "pH rule", "ODEs", "pH=7.4-a*Acid")

makeModel(glu)
simResults <- simulateModel(glu, times = seq(0,10,0.1))
plotResults(simResults, title="glu")

```

And the resulting plot is this where pH is increasing forever despite `Acid` is at a steady state concentration.

**Please note that the abovementioned code is just a simplification. It doesn't even look like aerobic glycolysis where 1 glucose produces lactate and 2 hydrogen atoms etc.**

I intend to include several acid- and base-creating reactions and evaluate the resulting pH. So could someone point me in the right direction about how to model the relationship between pH and acids or bases?

The answer does not need to be explicit code. I'll appreciate any solid comment on my design errors.

# Answer
> 2 votes
I figured that one has to include only the changes in the rules formula. Keeping the initial value of the `Species` or `Parameter` when defining circular or self-referencing formulas caused the error. So the correct rule becomes:

`pH=-a*Acid` instead of `pH=pH-a*Acid` or `pH=7.4-a*Acid`.

For convenience, I include the correct code crowded with several reactions producing both acids and bases, where pH gradually decreases as expected. The names are fictional, and just are better placeholders instead of x1 and x2, etc.

```
library("sysBio")

glu <- newModel("glu")
addSpecies(glu, "Glucose", 1000)
addSpecies(glu, "Lactate", 0)
addSpecies(glu, "Pyruvate", 0)
addSpecies(glu, "Oxalate", 0)
addSpecies(glu, "Acid", 0)
addSpecies(glu, "Base", 0)
addSpecies(glu, "pH", 7.4)

addMAreaction(glu, "Glucose = Lactate + Acid", "forwardL",
              "backwardL", name = "glycolysis")
addMAreaction(glu, "Lactate = Pyruvate + Acid", "forwardP",
              "backwardP", name = "lactolysis")
addMAreaction(glu, "Pyruvate = Oxalate + Base", "forwardO",
              "backwardO", name = "oxalate synthesis")

addMAreactRate(glu, "forwardL", "assigned",
               "k_glucose_to_lactate*Glucose")
addMAreactRate(glu, "backwardL", "assigned",
               "k_acid*pH+k_lactate_to_glucose*Lactate")
addMAreactRate(glu, "forwardP", "assigned",
               "k_lactate_to_pyruvate*Lactate")
addMAreactRate(glu, "backwardP", "assigned",
               "k_acid*pH+k_pyruvate_to_lactate*Pyruvate")
addMAreactRate(glu, "forwardO", "assigned",
               "k_pyruvate_to_oxalate*Pyruvate")
addMAreactRate(glu, "backwardO", "assigned",
               "k_base*pH+k_oxalate_to_pyruvate*Oxalate")

addParameters(glu, "k_glucose_to_lactate", 0.02)
addParameters(glu, "k_lactate_to_glucose", 0.001)
addParameters(glu, "k_lactate_to_pyruvate", 0.01)
addParameters(glu, "k_pyruvate_to_lactate", 0.0001)
addParameters(glu, "k_pyruvate_to_oxalate", 0.01)
addParameters(glu, "k_oxalate_to_pyruvate", 0.001)
addParameters(glu, "k_acid", 0.0001)
addParameters(glu, "k_base", 0.0001)

addRule(glu, "pH rule", "ODEs", "pH=-k_acid*Acid+k_base*Base")

makeModel(glu)
simResults <- simulateModel(glu, seq(0,3,0.01))
plotResults(simResults, title="glu")

```

The head of the resulting data.frame is:

```
  time   Glucose  Lactate   Pyruvate      Oxalate     Acid         Base       pH
1 0.00 1000.0000   0.0000  0.0000000 0.000000e+00   0.0000 0.000000e+00 7.400000
2 0.01  844.0001 155.0615  0.9384259 1.402282e-05 156.9384 1.402282e-05 7.399914
3 0.02  785.0869 210.3646  4.5476263 7.892761e-04 219.4615 7.892761e-04 7.399719
4 0.03  772.1557 218.6462  9.1924839 5.655403e-03 237.0425 5.655403e-03 7.399489
5 0.04  768.8212 217.2720 13.8876293 1.917807e-02 245.0856 1.917807e-02 7.399248
6 0.05  767.0383 214.4883 18.4278053 4.550630e-02 251.4350 4.550630e-02 7.399000

```

and resulting graph looks like:

---
Tags: r, modelling
---


# Question
Title: Get & Annotation of isomiR (isoform of miRNA) profile from RNA-seq
I should get and analyse isomiR (isoform of miRNA) profile from RNA-seq with FASTA format. The best tool for this is "DeAnnIso". But here is a problem with this tool, it only allows you to upload FASTA format files with a maximum size of 50MB. However, the size of my FASTA file is more than 1000MB. I tried to compress the file in different ways and uploading, but it doesn't work. Does anybody know any way to solving this problem?

Or Does anybody know **any other easy way to get isomiR profile from RNA-seq with FASTA format?**

A sample of RNA-seq like this:

```
>RR42357.95441.1 I-D00105:274:H7MDXX:1:112:8124:865 length=51
TGGAATTCTCGGGTGCCAAGGAACTCCAGTCACCTTGTAATCTCGTATGCC
>RR42357.95442.1 I-D00105:274:H7MDXX:1:112:8133:866 length=51
AAAGTGCTACTACTTTTGAGTCTTGGAATTCTCGGGTGCCAAGGAACTCCA
>RR42357.95443.1 I-D00105:274:H7MDXX:1:112:8144:8637 length=51

```

Thank you very very much!

# Answer
The obvious answer is to break up the file size and upload each separately, First I'd check the file size,

```
wc infile.fa

```

Output

> XXXXX YYYY ZZZZZ infile.fa

The number you want is XXXX - the number of lines in the file.

You'd then decide how many files you wanted to break it up into ... in your case divide XXXXX by 3. Then .. for example

```
sed -n '1,5000000p' infile.fa > outfile.fa
sed -n '5000001,2000000p' infile.fa > outfile2.fa
sed -n '2000001,4000000p' infile.fa > outfile3.fa

```

---

I'm not certain this would be cool however. An miRNA expert could answer this question very quickly. I think is you should check your data for sequence duplication. Simply piping the fasta file through a simple distance matrix and screening for "0" values would confirm this.

Creating a frequency table of duplicate RNA must be a feature of Bioconductor, but in any case screening for duplication (technically duplicate records) is a formal SQL database command.

It can be done a Perl/Python but its a bit clunky for a data set this size. I'd do it by setting the sequence as the key within a hash, the inefficiency has been solved (but its a weird way to do it).

You would have a database for "headers" for a given sequence, which would take up space, but I would assess the potential level of duplication first.
> 1 votes

# Answer
The tool you want to use accepts compressed files, so the first thing you should do is compress your file. I created a 111M file of random, 51bp sequences, each with a random, 52 character header and a random 51 nt sequence:

```
$ ls -lh file.fa
-rw-r--r-- 1 terdon terdon 111M Feb  9 12:21 file.fa

$ head file.fa
>YZF8SGLVBPCKMY18GJ94XTRCS61BB4OZKUZ61D 1J6R7H3 EJQA5 length=51
ACGAAGGGCATTTGGCAACGAACCTATACTATATCACGAGCAATTAGGCAG
>SPR0INGWKMILAXZZPF73BKRV73ELL1ROZSHYO1E9JX6P5X6MUZ20 length=51
TTAAGATTCCGAGAGAGTCTGTGGTTGGGTCCCTCTCCATAATGCTTTCAT
>U CB82L2M3BWJ9ZZZG3V6DL0Y0FTW8BXIFTU0TK7SL VN5QE42J6 length=51
TGGCTGGGAGGGATAGTAGCGAGATTATCCAGTCCGGTGTTCACAAGCCGG
>TGH62IM2TJK2QJTIVWGZR IVKOW94X7D52KGM7DJP61KKF81KPYA length=51
GCCGCAATACCAAAAAGATGGTGTGATTAGGAAATAAATAGGATCCAACAA
>L5VS9AVEF3UTEO3341STET092YEZ43ESVH32W7GI6Q27OIR NTXF length=51
CTATAGACGACTACGAACGACTGCAGGTGTGCCAGTTCTCGGTGAAACGGA

```

If we now compress that file:

```
$ tar cvzf file.fa.tar.gz file.fa
file.fa
$ ls -lh file.fa.tar.gz
-rw-r--r-- 1 terdon terdon 56M Feb  9 12:22 file.tar.gz

```

You will likely get a significantly better compression ratio since your headers won't be random and will likely be very repetitive.

If that isn't enough, you should try simplifying your headers. For instance, just use an increasing number:

```
$ perl -pe 's/^>.*/">" . ++$i/e' file.fa > file.short.fa
$ ls -lh file.short.fa
-rw-r--r-- 1 terdon terdon 58M Feb  9 12:23 file.short.fa

```

See what a huge difference that one change made? And if we now compress the file:

```
$ tar cvzf file.short.tar.gz file.short.fa
file.short.fa
$ ls -lh file.short.tar.gz
-rw-r--r-- 1 terdon terdon 19M Feb  9 12:24 file.short.tar.gz

```

Well below the threshold!

---

If compression isn't enough, you should make sure you have no duplicate sequences in the file. Use the `FastaToTbl` and `TblToFasta` scripts I have posted here and:

```
FastaToTbl file.fa | sort | uniq | TblToFasta > file.uniq.fa

```

Finally, if none of these helps enough, split your file into multiple smaller ones:

```
split -l 100000 file.fa

```

That will create multiple files of 100000 lines (so 50000 sequences, assuming all your sequences are one line) each.
> 1 votes

---
Tags: rna-seq, fasta
---


title: What is Feature Barcoding technology?

text:
# Question
Title: What is Feature Barcoding technology?
I heard the term Feature Barcoding a lot so I find resources to understand it. I have read docs of 10xgenomics and I see that `Feature Barcoding technology is a method for adding extra channels of information to cells by running single-cell gene expression in parallel with other assays`, this technology also do something with t-SNE and something with matrices.

I can't imagine the whole picture because I don't know what are those `adding extra channels of information to cells`, `single-cell gene expression` and `parallel with other assays` mean. And what is the input and output of this technology? Is there any paper go along with this technology?

# Answer
> 3 votes
"Feature Barcoding" refers to using antibodies linked to specific oligos. The antibodies can then bind to cell surface markers of interest and the oligos they're conjugated with turned into normal scRNA-seq reads. This effectively means you can get a FACs profile for each cell without needing to run it through a FACs sorter and individually label it.

The "feature barcoding" nomenclature comes from 10X I think, but there are a few similar protocols (maybe they're exactly the same in the end, I'm not sure about the details) such as CITE-seq and REAP-seq. You can see a talk from 10X on this here.

---
Tags: barcode, features
---


title: Installing DESeq2 in Ubuntu

text:
# Question
Title: Installing DESeq2 in Ubuntu
I am trying to install DESeq2 in my Ubuntu with R version 3.5.1. According to the package repository in Bioconductor the version should be 3.5.

```
> R.version

platform       x86_64-pc-linux-gnu
arch           x86_64
os             linux-gnu
system         x86_64, linux-gnu
status
major          3
minor          5.1
year           2018
month          07
day            02
svn rev        74947
language       R
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

```

then I attempt to install DESeq2, following the instructions in Bioconductor:

```
> if (!requireNamespace("BiocManager", quietly = TRUE))
       install.packages("BiocManager")
   BiocManager::install("DESeq2", version = "3.8"

```

and at the end I get:

```
Warning messages:
1: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘RCurl’ had non-zero exit status
2: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomeInfoDb’ had non-zero exit status
3: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘annotate’ had non-zero exit status
4: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘GenomicRanges’ had non-zero exit status
5: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘genefilter’ had non-zero exit status
6: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘geneplotter’ had non-zero exit status
7: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘SummarizedExperiment’ had non-zero exit status
8: In install.packages(pkgs = doing, lib = lib, repos = repos, ...) :
  installation of package ‘DESeq2’ had non-zero exit status

```

I have visited $2\*10^{9999}$ discussions on this issue and none of them helped me solving the problem. In addition, I didn't have any problem installing other packages like limma and edgeR, so I assume that this might be specific to DESeq2 or one of it's dependencies.

Thank you for the help!

# Answer
First I installed this libraries:

```
sudo apt-get install libcurl4-openssl-dev libxml2-dev

```

I then installed RCurl:

```
BiocManager::install("RCurl")

```

and lastly:

```
BiocManager::install("DESeq2", version = "3.8")

```
> 4 votes

# Answer
Have you tried just typing,

```
install.packages("BiocManager")

```

Without all the proceeding if statement? If this doesn't install correctly nothing else will. Personally I would be tempted install each of the named libraries independently. If that works okay, e.g. its not a problem with the mirrors then I'd contact the developer michaelisaiahlove at gmail.com
> 1 votes

---
Tags: deseq2, software-installation, linux
---


# Question
Title: Question about Co-expression analysis and finding targets for lncRNAs
I have a dataset with 88 tumor and 113 normal samples. Among the 50k genes after filtering there are a total of 28k genes (both mRNAs and lncRNAs)

I wanted to do co-expression analysis between lncRNAs and mRNAs to find the targets for each lncRNA. I'm thinking to use WGCNA for co-expression network analysis.

Few doubts:

1) For Co-expression network analysis (WGCNA) should I use only tumor samples (or) both tumor and normal samples?

2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?

3) Heard that with bedtools we can know which genes are 10kb upstream/ downstream of lncRNAs to find neighboring pc genes. But no idea how to do that. Could you please give me an example for this.

Thank you

# Answer
You have several options to approach this with WGCNA (weighted correlation network analysis). You can run a WGCNA on the combined set, identify modules and select those lncRNAs for further follow-up that have relatively high intra-modular connectivity (i.e., are intramodular hub genes). I would suggest this approach first.

Second option is to run WGCNA on mRNA for coding genes only and screen module eigengenes and lncRNAs for association with one another. This approach is usually helpful when the two data types (mRNA and lncRNAs) are different which I suppose is not really your case (mRNA and lncRNAs were probably measured together and are overall similar).

I would start using all samples together (tumor and normal). Hopefully you'll find some interesting and plausible information; in the second pass I'd perhaps create separate (or consensus) networks for tumor and normal samples.

I am not familiar with bedtools, but you should have some sort of gene/transcript annotation with your data (it was surely aligned to a reference genome and summarized to some sort of annotated set of genes). The annotation should contain chromosome and base pair position, which would make it rather easy to find genes within a specific interval.
> 4 votes

# Answer
A cursory search to bedtools documentation will reveal the `bedtools closest` feature - which might be exactly what you are looking in your third question. You can download the list of protein-coding genes from Ensembl biomart or UCSC table browser in bed format, and also convert the list of your lncRNAs into a bed compatible file. Rest is easy:

```
bedtools closest -a lncRNAs.bed -b proteinCoding.bed [-iu/-id] -D ref

```
> 2 votes

# Answer
> 2) If not WGCNA, if I use cor function in R with method Pearson, on what cutoff I should select the target genes?

Good call, this is a nice approach. It is a good question because the central issue is that a purists approach of 0.7 results in insufficient information to understand the behaviour of the data.

The guys here were using 0.4, in my experience co-variance of 0.5 and above is good evidence, whilst 0.7+ is strong evidence.

I've done this (not on RNAseq) and I used 0.25 to 0.45 and 0.45 and above then used a scatterplot with one half of the graph showing low co-variance and the other half strong co-variance. So in one sense I ignored my own advice (0.5+). In my defence its easy to reassign the boundaries. Each scatterplot (1000s) was then subject to tSNE analysis to assess the group behaviour. Its a different context but gives an idea.

Networks are cool for visualisation but they are not particularly quantitive for most (not all) network theory.

---

Just to explain to everyone else, what happens if you get boundary effects with the data and you are assigning a discrete boundary on a continous distribution (its a negative binomial distribution in my case). I observed a vast amount of results occurring between 0.2 to 0.25 - but 0.2 isn't good evidence, so that data has to be lost. These sort of boundaries were in days before genomic throughput. Thus in my context I didn't consider a hard boundary particularly informative.
> 1 votes

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


# Question
Title: What is the index fastq file (sample_I*.fastq.gz) generated when demultiplexing Illumina paired-end runs?
What is the index fastq file that comes with some Illumina sequencing datasets? (The `samplename_I*.fastq.gz` file.)

For example, I recently received some 10X Chromium reads for two libraries sequenced on the same lane. This was a 2x150 sequencing run, so there should be two fastq files. The sequencing center demultiplexed the libraries and generated two separate directories - one for each library.

* In each directory there are three fastq files:
    + `Mysample_I1_001.fastq.gz`
    + `Mysample_R1_001.fastq.gz`
    + `Mysample_R2_001.fastq.gz`

I know that the R1 fastq file contains the forward reads and the R2 fastq file contains the reverse reads, but what is the index file? The first few sequences do not match the 10X indices for this library (`ACATTACT`, `CAGCCCAC`, `GGCAATGG`, `TTTGGGTA`).

```
>bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | head
NTTGGGTA
NGCAATGG
NAGCCCAC
NAGCCCAC
NCATTACT
NCATTACT
NGCAATGG
NAGCCCAC
NTTGGGTA
NAGCCCAC

```

# Answer
# tldr - The I\*.fastq.gz file contains the read index sequences.

## long explanation

Illumina uses a program called `bcl2fastq` to demultiplex sequencing runs.

This software takes a list of samples and their associated indices and uses those sequences to make one or more fastq files per sample, binned by one or two index sequences on either end of the sequencing molecule (i5 and i7 indices, see page 6 for HiSeq).

Illumina sequencing isn't perfect though, and there are sometimes errors in reading the index sequence. For example this index, `CAGCCCAC` can easily have read errors in the `A` sandwiched between many `C`s: `CAGCCCAC -> CAGCCCCC`.

Instead of throwing out all of the reads with indexing sequencing errors, the `bcl2fastq` program includes reads that are discernibly derived from a true index for that sample as long as there is no overlap with another sample.

In the example provided above, the four indices for `mysample` were: `ACATTACT`, `CAGCCCAC`, `GGCAATGG`, `TTTGGGTA`. If we look at all of the index sequences in the `Mysample_I1_001.fastq.gz` file, we will see that sequences with the four correct indices are indeed the most prevalent, however there are also reads that have indices with sequencing errors derived from the correct indices.

```
bioawk -cfastx '{print($seq)}' Mysample_I1_001.fastq.gz | \
    sort | uniq -c | sort -k1 -nr | head

count      index       source
41362311   CAGCCCAC    True index
37209190   GGCAATGG    True index
36863213   ACATTACT    True index
33674467   TTTGGGTA    True index
 1140358   NAGCCCAC    CAGCCCAC
 1026099   NGCAATGG    GGCAATGG
 1016754   NCATTACT    ACATTACT
  933342   NTTGGGTA    TTTGGGTA
  119626   TTTGGGGA    TTTGGGTA
   98657   GTTGGGTA    TTTGGGTA
   96625   GGCAATGA    GGCAATGG

```
> 9 votes

# Answer
It took me a while to figure out that the "**index**" is the same thing as the "**barcode**" that says which sample each sequence is from on a multiplexed run.

If your data is *not demultiplexed* (single `R1.fastq` and `R2.fastq` files contain the information for multiple samples), then this `I1.fastq` file is what you use to assign each sequence to a sample (ie to "demultiplex"). The spreadsheet listing which barcode/index is which sample is often called a "mapping file."

If your data is *already demultiplexed* (you have separate `R1.fastq` and `R2.fastq` files for each sample), then you don't really need it... however some people use the index sequences in their quality control workflows. For example, this 2016 paper in BMC Genomics by Wright and Vetsigien "Quality filtering of Illumina index reads mitigates sample cross-talk": https://doi.org/10.1186/s12864-016-3217-x
> 2 votes

# Answer
I looked at 10x data a couple of years ago, and at that time, their protocol would use four separate indexes, with balanced index sequences, as a way of introducing diversity into the library. These four indexes were separate from the individual oligos used to mark single cells. For an individual sample, you would need to combine the reads from all four separate indexes, which is how it looks like in your data. There is a brief explanation here: https://support.10xgenomics.com/single-cell-gene-expression/index/doc/specifications-sample-index-sets-for-single-cell-3
> 0 votes

---
Tags: ngs, illumina, 10x-genomics, demultiplexing
---


# Question
Title: Error in as.vector(x) : no method for coercing this S4 class to a vector
I tried to run the following code in R studio. Everything worked fine, except at the last step \[write.table(mdat, "recount\_mdat.csv")\] when I tried to export the 'mdat', I got the following error:

```
Error in as.vector(x) : no method for coercing this S4 class to a vector

```

My code:

```
 library('dplyr')
 library('recount')
 library('magrittr')
 library('ffpe')
 library('RSkittleBrewer')
 library('SummarizedExperiment')
 library('devtools')
 trop <- RSkittleBrewer::RSkittleBrewer('tropical')
 colon_proj <- c('SRP025982', 'SRP42161')
 if(any(!file.exists(file.path(colon_proj, 'rse_gene.Rdata')))) {
 sapply(colon_proj, download_study)
 }

 blood_proj <- c('SRP066834', 'SRP041736', 'SRP060416')
 if(any(!file.exists(file.path(blood_proj, 'rse_gene.Rdata')))) {
 sapply(blood_proj, download_study)
 }
 proj <- c(colon_proj,blood_proj)

 dat <- lapply(proj, function(x) {
 load(file.path(x, 'rse_gene.Rdata'))
 return(rse_gene)
 })
 proj

 sapply(dat, dim)

 metadata <- all_metadata('sra')

 write.table(metadata, "recount_metadata.csv")

 mdat <- do.call(cbind, dat)

 write.table(mdat, "recount_mdat.csv")

is(mdat)
## "RangedSummarizedExperiment" "SummarizedExperiment" "Vector" "Annotated"

```

# Answer
You can't cbind a bunch of obscure object types.

If you want merged count tables you should do this:

```
mdat <- do.call(cbind,lapply(dat,assay))

```

Where **row.names** are *Ensembl gene IDs* and **col.names** are the *SRR accessions*.

Then run your table writing command.

-

If you want the coordinates of your genes then do this to make a bed with the genomic locations:

```
library(rtracklayer)
export.bed(rowRanges(dat[[1]]),'gene_coords.bed')

```
> 6 votes

---
Tags: rna-seq, bioconductor, rstudio
---


# Question
Title: Using preprocessing/alignment functions on the server
I am new to bash and the processes behind cluster computing in general and need some help with understanding some basics.

After looking all over the internet and this forum (+ askUbuntu) I found nothing that addresses this issue.

I have a series of raw RNA sequencing files: .bam and wish to utilize samtools (or others) to begin the data preprocessing steps defined in many workflows.

I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error: `-bash-samtools: command not found`

When I try to install the package, say samtools, by running `sudo apt install samtools`, I get an error stating I do not have permission to do such.

What do I do? Must I bring up this issue to the server provider or is there a way around it?

Thank you for your time,

# Answer
> 3 votes
The short answer is to use conda. In the bioconda channel we have most of the tools used in bioinformatics, such as samtools. You do not need administrator permissions to install conda or packages with it, so your lack of `sudo` ability is not an issue.

> I have installed samtools and many other programs necessary for the workflow, however, once I am on the remote server in bash, and try and run a function over the files of interest, it gives me the error

Installing packages on your local machine won't affect the remote server unless they share a file system, so this is expected.

As an aside, BAM files aren't really "raw", they're considered processed data (fastq/fast5 files are considered truly "raw" data in the sequencing world).

# Answer
> 2 votes
At a very wild guess

```
-bash-samtools

```

This looks shaky. For a start you need a space between the command and the application in all Linux/Unix/Mac OS X applications (its complicated to explain right now \[see below\] just accept it).

---

Wait you are using `sudo` and thats how you are installing on a local machine. `sudo` on a remote machine will never work, because you are not the systems admin. `sudo` on a local machine is not cool (see below)

---

Anyway I need to do some next-gen and have 6 MiSeq kits waiting for action so I installed it from here: http://www.htslib.org/download/

Okay so here's my tcsh history

```
 4  21:44   cd sam* # like I was in the Desktop
 7  21:44   more README # I needed to know how they wanted the installation
 8  21:45   ./configure # okay there are setting the paths for me - saves me a job
11  21:46   make # like fairytale stuff for a load of C scripts
16  21:49   ls | grep samtools # I wanna know if its compiled and I see samtools* ... yep
17  21:49   ln -s ./samtools* ~/bin # I'm gonna pop a link in my bin so after I type
18  21:49   rehash # 'cause I still insist on tcsh rather than bash (don't do it in bash)

```

However, the following is a bit creepy because I get loads of perl scripts in my /usr/local/bin ...

```
make install

```

Anyway,

```
samtools

```

It took 4 mins on an OS X machine, using tsch, after I just quaffed 2 pints at my pub (10pm London time). So probability is set against me ....

**OUTPUT**

```
Program: samtools (Tools for alignments in the SAM format)
Version: 0.1.18 (r982:295)

Usage:   samtools <command> [options]

Command: view        SAM<->BAM conversion
         sort        sort alignment file
         mpileup     multi-way pileup
         depth       compute the depth
         faidx       index/extract FASTA
         tview       text alignment viewer
         index       index alignment
         idxstats    BAM index stats (r595 or later)

```

etc ... BTW seriously, fun over, avoid "`sudo`" if you can, for lots of reasons (which I'm not gonna explain). Its a last resort and usually means your paths are incorrectly set.

sudo will be refused on a remote machine, without question. You COULD ask systems admin to do this for you. HOWEVER, I'd try and install yourself without `sudo` and your bioinformatics will benefit.

---

```
Program: samtools (Tools for alignments in the SAM format)
Version: 1.9 (using htslib 1.9)

Usage:   samtools <command> [options]

Commands:
  -- Indexing
     dict           create a sequence dictionary file
     faidx          index/extract FASTA
     fqidx          index/extract FASTQ
     index          index alignment

```

**Output** following `make install`

---
Tags: rna-seq, scrnaseq, samtools, bash
---


# Question
Title: Error in PSI-BLAST when running psipred locally
I am very new to bioinformatics and Linux and am trying to run psipred locally on Ubuntu 18.04 but cannot get PSI-BLAST to work. I have downloaded BLAST legacy version 2.2.26 and unpacked it, uniref90 database and psipred version 4.01

I installed psipred according to instructions:

```
cd psipred/src
make
make install

```

And this seems to work fine

I used the following code to create a database, which proceeded without errors:

```
makeblastdb -in uniref90.fasta -dbtype prot -out uniref90

```

In runpsipred I changed the directories:

```
# The name of the BLAST data bank
set dbname = /home/hilda/uniref90

# Where the NCBI programs have been installed
set ncbidir = /home/hilda/blast/bin

# Where the PSIPRED V4 programs have been installed
set execdir = ./bin

# Where the PSIPRED V4 data files have been installed
set datadir = ./data

```

I then try to run psipred:

```
tsch ./runpsipred example/example.fasta

```

And get the error

```
Running PSI-BLAST with sequence example/example.fasta ...
FATAL: Error whilst running blastpgp - script terminated!

```

I am not sure how to solve this as I do not understand what is going wrong. As I said I am new to Linux. I use the command tsch to run psipred (otherwise I get the error "command not found"), but not to install the database, could this be a problem? Or have I missed something in the setup process? Running psipred does create a copy of the example.fasta in the psipred folder named psitmp4286007f010.fasta, but that is all that seems to happen.

I would greatly appreciate any suggestions on how to make it work!

Edit: I tried to run with BLAST+ (which seemed to work without tsch) and get the following error:

```
BLAST Database error: No alias or index file found for protein database [/home/hilda/uniref90] in search path [/home/hilda/psipred::]

```

I have also tried setting dbname to home/hilda/uniref90/uniref90 and got the same results. there is a uniref90.pal file in the uniref90 directory, which I understand is the alias file? Other forum answers seem to use pfilt to index (?) the database, but there is no pfilt in my psipred directory. Any ideas how to fix the problem?

# Answer
This error could be due to the fact that you are using legacy blast, have you tried using BLAST+ instead? If you wish to use PSIPRED with BLAST+, then use the runpsipredplus script here, rather than the normal runpsipred script.
> 4 votes

---
Tags: software-installation, psi-blast, makeblastdb
---


# Question
Title: How to quantify similarity of genomes and find differences in set of S aureus genomes?
I have around 500 annotated proteomes of different bacterial strains and would like to quantify their similarity (or difference). I found `gt genomediff` from genometools gives me some scores that I can use to generate nice clusters, but I am not sure whether that tools really works. The fasta-files that I use contain multiple sequences.

I ran some tests and it looks ok.

```
1_reference.fna:
>1
TAAGTTACT
>2
TAAGTTACA

2_eq_to_ref.fna:
>1
TAAGTTACT
>2
TAAGTTACA

3_tags_diff.fna:
>1asdfadf
TAAGTTACT
>2asdfasdfffa
TAAGTTACA

4_orde_diff.fna:
>2
TAAGTTACA
>1
TAAGTTACT

5_add_subse.fna:
>1
TAAGTTACT
>2
TAAGTTACA
>2
TTACA

6_point_mut.fna:
>1
AAAGTTACT
>2
TAAGTTACA

7_different.fna:
>1
TAAGTTACT
ATTACCTAA
>2
AAAAAAAAA

```

Then:

```
gt genomediff --indexname test *fna
7
1_reference.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
2_eq_to_ref.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
3_tags_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
4_orde_diff.fna 0.000000        0.000000        0.000000        0.000000        0.206969        0.199527        0.794496
5_add_subse.fna 0.206969        0.206969        0.206969        0.206969        0.000000        0.212596        1.349206
6_point_mut.fna 0.199527        0.199527        0.199527        0.199527        0.212596        0.000000        0.569180
7_different.fna 0.794496        0.794496        0.794496        0.794496        1.349206        0.569180        0.000000

```

I am running the clustering on the difference matrix calculated with genomediff: "These distances are Jukes-Cantor corrected divergence between the pairs of genomes, that is, the number of mutations per base between them."

Currently, we are studying S aureus. We genomes are assembled genomes (three different methods). My guess is that the sequences from the plasmids are present. Furthermore, we do have drug-resistance measured in culture. So, we will be able to compare the genomes and the resistances.

# Answer
... stuff of original post deleted. On second thoughts what you might be doing is a template assembly of your genomes. It is a possible interpretation of the 1.2. fasta sequences above (i.e. 1 is template).

A microbial genomics professional would advise also performing a de novo assembly, particularly if you are interested in presence and absence. The reason is that if a gene is present in your query and absent in your template, you will miss it. Again it all depends on what bacteria are being assessed, some are more prone to "genetic islands" than others.

You need a collaborator beyond that.

For the onlookers here bacterial genetic behaviour is very different to eukaryotes and what they get up, switching DNA etc.., to would appear bizzare from a higher eukaryotic world.

---

You mentioned a comparative analysis of 500 strains: I have worked on this in MRSA. Note, the bacteria is important considering the approach you will adopt.

Anyway, you want a single aligned file to produce a phylogeny, Bayes or likelihood. This is a model of point mutations. Bootstrapping doesn't really help because of SNP differences between isolates.

There is a complex problem. Generic phylogeny can hit the buffers in my opinion because often:

* low numbers of SNP differences between isolates.
* The error across the genome is poorly defined
* Multi-clonal infection ignored and your isolates not cloned.
* The other issue is the QC of the genome, I've been fairly astonished by the variation (it affects the tree)

Put it all together and you can get a nice tree but the topology of a given MLST is in my opinion unlikely to be correct. A referee may not bother with this, but they might.

You then map the phenotype against the phylogeny (the tips of the tree) and look for clusters. HOWEVER, there are problems:

* Drug-resistant e.g. mec cassette (methicillin-resistant ) are often on plasmids and these are lost during isolation. Albeit they can integrate onto the genome. So an isolate can be drug resistance but you fail to find the gene.
* Drug-resistance transmission is unlikely to cluster, its too quick, so you don't see nice tight clusters against a theoretically perfect tree
* The best approach is drug-resistance in culture.

Hierarchical clustering, has many meanings. It is used in presence/absence data, I demonstrated the method failed (on bacteria) based on cluster analysis. It is currently revised and used in unsupervised deep (or machine) learning prior a training method. The question then is what is what are you modelling? Presence/absence (of genes), point mutations, epidemiological data, drug-resistance/non-drug resistance?

My assessment is that given the potential complexity of these biological scenarios and the quality of the question you should seek formal collaboration and this is before considering what bacteria are being assessed. Some of epidemic drug-resistent bacteria rip up this rule book.
> 4 votes

---
Tags: phylogenetics, sequence-alignment, phylogeny, bacteria, genometools
---


# Question
Title: Questions regarding Nanopore sequencing analysis
I am new to Nanopore sequencing analysis. I have a couple of questions regarding it which are as follows:

1. How do I know if my fast5 file is multiread or single read file?
2. Is there a way to combine all the fast5 files into a single fast5 file?
3. How can I extract the flowcell type and sequencing kit information from the fast5 file? Is there a command that I can use to extract it? I know I can use the following command to retrieve it:

```
h5dump -a /UniqueGlobalKey/context_tags/sequencing_kit -a /UniqueGlobalKey/context_tags/flowcell_type <fast5 file>

```

However, it works well for single read file and not for multiread files. Please correct me if I am wrong. Thanks!

# Answer
*\[Note: As of 2024, ONT now uses pod5 as the standard output format from their MinKNOW sequencing program, and *all* pod5 files are "multi" pod5 files.\]*

**1.** The quickest way to know if it's a single or multi fast5 is the file name. If it mentions channels and read numbers (e.g. `..._read_9_ch_471_...`, then it's a single read file. Alternatively, the multi-read files have only read IDs at the top level (i.e. read\_\<32-character-string-with-additional-dashes\>). Example:

```
$ h5ls <read_file_name>.fast5 | head
read_0115573a-042a-42a3-b4e2-f0338b4a6c66 Group
read_015b870d-8a03-4f22-a73e-fb9f833b94c7 Group
read_0167393a-dede-4a27-91af-494622b63261 Group
read_01e658de-853b-4f58-93ad-87fbd735254c Group
read_0202acc2-33f4-4d29-9dac-c7f80c758c86 Group
read_022d68f0-6f8e-4124-bfb0-11262efa4a67 Group
read_0258a3c9-4ffe-42ab-a42d-89466fcefbf5 Group
read_032da6db-e84a-43ec-8c25-bcaead992066 Group
read_03383754-be43-420b-87fc-60c2eb833963 Group
read_036b4824-bb0a-42dd-ade0-29d89112ff07 Group

```

**2.** ONT provide a `single_to_multi` program for converting to the multi-read format. Details can be found on this github page. Here's a quick example of how it could be run:

```
$ single_to_multi_fast5 --input_path fast5 --save_path fast5_multi \
    --filename_base multi_file_myRunID --batch_size 4000 --recursive

```

**3a.** For multi-fast5 files, the kit/cell attributes are in the `context_tags` group of the file (because a multi-fast5 file could contain fast5 elements from multiple different sequencing runs with different run parameters):

```
$ h5dump -g /read_<readID>/context_tags multi_file_myRunID_0.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
   ATTRIBUTE "flowcell_type" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "flo-min106"
      }
   }
--
   ATTRIBUTE "sequencing_kit" {
      DATATYPE  H5T_STRING {
         STRSIZE 11;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "sqk-pbk004"
      }
   }

```

**3b.** For single fast5 files (at least, in the most recent version of them), they are found in the `UniqueGlobalKey` group:

```
$ h5dump -g /UniqueGlobalKey single_file_read_X_ch_Y_strand.fast5 | \
  grep -A 11 -e 'sequencing_kit' -e 'flowcell_type'
      ATTRIBUTE "flowcell_type" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "flo-min106"
         }
      }
--
      ATTRIBUTE "sequencing_kit" {
         DATATYPE  H5T_STRING {
            STRSIZE 11;
            STRPAD H5T_STR_NULLTERM;
            CSET H5T_CSET_ASCII;
            CTYPE H5T_C_S1;
         }
         DATASPACE  SCALAR
         DATA {
         (0): "sqk-pcs108"
         }
      }

```

\[note that this file structure may change in the future\]
> 6 votes

---
Tags: nanopore, long-reads
---


# Question
Title: Targeted NGS, up to 99% of reads have been marked as duplicates
Currently I'm performing whole analysis (pipeline from \*.fastq to \*.vcf) of 41 samples (targeted NGS). I rely on GATK best practices, however with some modifications. I decided to use the following tools:

```
#mapping
bwa mem (with mem - alternate haplotypes)

#preprocessing
RevertSam (Picard)
MergeBamAlignment (Picard)
MarkDuplicates (Picard)
BaseRecalibrator (GATK)
ApplyBQSR (GATK)
AnalyzeCovariates (GATK)
SortSam (Picard)

#call, filter and annotate variants
HaplotypeCaller (GATK)
GenotypeGVCFs (GATK)
snpEff, SnpSift

```

Nevertheless, after marking the duplicates with `MarkDuplicates` about 99% of reads have been marked as duplicated, hence the coverage decreased from several hundreds to literally up to five. So I decided to omit this step.

Have you ever encountered a similar problem? What is your opinion about this?

Best Regards, Adam

# EDIT

The Agilent Haloplex custom target enrichment kit (HaloPlex Target Enrichment for Illumina Kit) has been used to prepare the libraries, sequencing done using Illumina HiSeq 2500. Targeted NGS of all exons (only) of 90 genes.

Here I provide how the final bam look like. This is the result without duplicate marking, otherwise there're only several reads per exon...

# Answer
> 1 votes
Most duplicate marking programs do it by looking at the start and end location of hits. In this case, if two reads have identical start and end locations for both reads, then they will be considered \[PCR\] duplicates of each other. When running through a duplicate filter, you will only get one read for each unique pair.

This is much more likely for targeted sequencing, because the length covered by the template is much smaller. It's also more likely that you actually *do* have lots of PCR duplicates for the same reason.

# Answer
> 3 votes
If you're doing target capture, especially of a small region, it is entirely possible that there are only five unique reads. Target capture is molecularly finicky and sometimes the entire region is not captured.

---
Tags: gatk, markduplicate
---


# Question
Title: How can I get plant Genome size, total number of genes, EST, and chromosomes number
How can I get plant Genome size, total number of genes, EST and chromosome number from phytozom and plantDB?

# Answer
> 2 votes
Ensembl plants has a REST API, you might be able to get all these details (except possibly EST numbers) there: http://ensemblgenomes.org/info/access/rest

For example this query gives you the Arabidopsis genome size: http://rest.ensemblgenomes.org/info/genomes/arabidopsis\_thaliana?content-type=application/json

---
Tags: public-databases
---


# Question
Title: Retrieve RNA sequencing data for human p53 colon cancer cell lines
I want to download the next generation RNA sequencing data in different cell lines. I want the data (alongwith normal) for p53 wild type, p53 knock-out and p53 null cell lines (human colon cancer). From where do I download this data? In fastq format if possible.

Currently, I am going through SRA, to look for the datasets, using their search engine. But the process is extremely tedious.

Also can someone please hint about how to download the data from SRA, in case that's the only way to get data. I'm getting very confused.

# Answer
I suggest you to try difference resources, including ArrayExpress. A simple search at ArrayExpress returns multiple experiments which might be of your interest. Just a pointer, ArrayExpress is also synced with GEO.
> 2 votes

# Answer
I agree with Kamil on the downloading method. Moreover, SRA will give you your file in the raw format which is usually a **Binary Alignment Map (BAM)** or **Sequence Alignment Map (SAM)** Download SRA toolkit available online, to convert your raw reads to fastq format. Follow the steps using command line. It is a very simple command to convert your SRA file to fastq which uses the *fastq-dump* command.

Further reading: https://www.ncbi.nlm.nih.gov/books/NBK158900/

Hope it helps!
> 1 votes

---
Tags: sequence-analysis, rna
---


title: Randomness in BLAST

text:
# Question
Title: Randomness in BLAST
So if you see the BLAST parameters it says

> The Expected value E is a parameter that describes the number of hits one can "expect" to get by chance when searching a database of particular size. It decreases exponentially as the score (S) increases. Essential, the E value describes the randomness background noise. For example an E value of 1 assigned to a hit can be interpreted as meaning that a database of the current size one might expect to see 1 match with a similar score simply by chance. Lower the E value the better."

How could there be something "BY CHANCE" it's a software right ?

# Answer
> 7 votes
The Blast E is the expected frequency of obtaining false positive. This will depend on your query size (number of nucleotide or amino acid residues) and the size of the database. With a short query and a large database you are more likely to have a sequence in your database that matches your query by simple chance. In summary, the smaller the query and the larger the database the greater the chance of a spurious result and that is what E is measuring.

# Answer
> 5 votes
Sure, so you start out with what is called a bitscore, which is a normalized to the score calculated from the alignment between the 2 seqs which depends on the following equation. It is **independent** of database size:

> (lambda * S - ln(k))/(ln)2

Then, the p-value of a local blast is just:

> 1/2^bitscore

So if your bitscore is 15, you need 1/32768 alignments before you will get a score as good or better (highly similar sequences) **BY CHANCE ALONE**. This addresses your "BY CHANCE" question.

Earlier I said the bitscore was **independent** of database size. The E value is just the p-value above normalized to the database size (so it is **dependent** on the db size) by the following equation:

> (query length * database length * p-value)

which simplifies a bit to:

> E = (query length\* (db length/2^bitscore)

I hope this helps!

# Answer
> 0 votes
Calculation of E for Blast is a good question.

Two methods:

* Use a Poison or binomial distribution.
* Use randomization

I think its no. 1 because no. 2 would use too much computational power. I did know but someone borrowed by O'Reilly book on the subject.

One of the contributors here will know, because I suspect their genomics developers.

---
Tags: blast
---


# Question
Title: Finding the members in a confusion matrix
I have inferred a confusion matrix of training and test set by neural network. I want to know which members are in the confusion matrix.

```
> table(test$TRG,predicted.nn.values$net.result)

    0 1
  0 5 4
  1 5 3
>

> head(train[,1:4])
  TRG      ALB     AQP9   CALML5
1   1 5.865827 8.190945 6.705303
2   1 6.998435 8.424261 8.505591
3   1 7.424512 8.716471 7.249556
4   1 7.442049 8.325263 8.809286
6   1 5.893411 8.199990 6.677618
7   1 7.288030 9.143510 7.088598
>

> head(test[,1:4])
   TRG      ALB     AQP9    CALML5
5    1 6.369307 7.954310  6.920290
8    1 6.181902 8.651442  7.225389
10   1 6.119359 9.345270  6.829623
13   1 7.775533 9.016272  7.976813
14   1 5.913656 9.484457 10.609013
18   1 6.603138 7.908560  8.827173
>

```

How can I know which patients are false positive or true negative? For example I have 5 patients predicted correctly, how can I know the name of these patients?

# Answer
This does not look a great confusion matrix, but it is a very, very cool approach. The cool thing about machine learning is the ability to examine lots of models within confusion matrices.

---

To answer your question, you want a heat map of the regression weights to establish who is who, i.e. which patient is false positive etc... This is the only way I'm aware of to analyse a neural network. This looks like R, I only know Scikit learn.
> 1 votes

---
Tags: r, machine-learning
---


# Question
Title: How to cluster the human genes by pathways/system-biology/metabolic properties?
I would like to make a chord plot from my data. I have a list of genes that according to my experiments are divided into 64 clusters of enrichment pattern.

I would combine my 64 clusters with a second clustering. A clustering of all the human genes into subgroups of "interacting" genes (e.g. both are pathway, co-expressed) will allow me to combine the two clustering data on a chord plot.

In the analogy of the following graph the continents will be my 64 clusters of the genes and the edges/chords will be assigned based on the biological function clustering. I guess that such a functional clustering can be done ore has been done using platforms such as: KEGG, genescf, or reactome.org.

1. Is there a biological pathway clustering of genes for lymphoblastoid or breast cell lines?
2. Are there any suggestions how to build the clustering between genes whose functions are strongly correlated in the cell?

**Edit**: the question is about getting the clustering from a system-biological perspective not how to plot chord diagrams.

Is there such a metric to quantify how likely two genes are correlated in their enrichment, function etc. For example, using STRING we can see that PIK3CA and PTEN are more co-functioning than PIK3CA and SF3B1. How can I do something like String on genome wide scale.

# Answer
> 2 votes
To measure if two genes are functional similar I developed the BioCor package. It calculates a similarity score between genes by the measuring the amount of pathways shared. However, it doesn't take into account the strength of the connection between two proteins. You could multiply the similarity by a value between 0 and 1 according to the number of evidences of the relation between the two proteins.

It is optimized for large amounts of data, you could even perform them in parallel, but I doubt you find information about pathways for most of the genes (In my informal test, I think that I could calculate similarities for about 1600 genes not the ~20000 of the human genome).

When I am interested in co-expression data I calculate the co-expression by correlation and add the similarity of the genes calculated with BioCor. The resulting distance matrix can be used for clustering. Maybe the section about using weighted co-expression network analysis and BioCor can be helpful.

To have a "biological pathway clustering of genes for lymphoblastoid or breast cell lines", you would need to have a specific signature/pathway of the cell lines. Currently I am not aware of accepted standard signature database. The single cell atlas could be useful to find expression of your cells, but I don't know if they provide a signature for cell line

# Answer
> 1 votes
Production of a chord diagram is described in detail here. You use the R package with either,

```
library(chorddiag)
library (circulize)

```

The first level clusters are predefined and the second layer of clustering is what you need to work out. The obvious stat is correlation, secondly (better approach IMO) you might look at covariance, although I don't exactly know what this "interaction" constitutes, which makes it difficult.

My personal view is the diagramatic representation doesn't really matter, its the robustness of the stats behind the diagram that counts.

I was a bit surprised to learn ggplot2 doesn't do circular diagrams :-(

> Lymphoblastoid Cell lines

I think you need to clarify the experiment on these cells lines, I assume its a time course, i.e. time series.

---

At a guess this is RNA expression, if you are using a chip then its between predetermined genes (i.e. spotted on the chip). If it is RNA seq then it can be anything. I suspect you are treating the cells and monitoring changes in RNA expression as a function of their biological (immunological?) activity. If the basic measure is time, i.e. a time series - thus two genes being expressed in the same unit time the best opening statistic is covariance. If all you want to do is make a chord diagram out of it, I think thats fine, this will form a 2x2 variance matrix for which any standard clustering will work (Jaccard distances \[there will be better ones\]) or direct plotting from the matrix. If you looking at unique expression behaviour to attempt to unlock some sort of pathway then you are into more serious regression time series (autoregression that sort of thing - this is complicated). If its oscillations then its a separate branch of stats based on sinosoids (sine waves) called periodicity analysis. HOWEVER, if the interaction is genetic similarity (gene families) e.g. mRNA then you want a phylogenetics measure, which if its clustering neighbor-joining is fine.

---

My personal view if it is the experiment I think, I would,

1. Construct a matrix of the covariance of the (I assume) time series
2. Make sure the covariance is standarised between -1 to 1 (Pearson approach), R-cran.
3. Assign a tight cut-off (tricky because there's no "hard cut-off") but say 0.7 (you need to experiment a bit)
4. Plot the raw covariance on the chord diagram.

Easy, simple and highly visual.

---
Tags: ngs, gene-expression, pathway, human, genetics
---


title: Determining Read Groups

text:
# Question
Title: Determining Read Groups
Which Read Groups are correct:

```
java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84

```

or

```
java -Xmx4G -jar picard.jar AddOrReplaceReadGroups \
I=$SNIC_TMP/WCRO84_S23_L005.bam \
O=WCRO84_S23_L005.bam \
RGID=WCRO84_S23_L005 \
RGLB=WCRO84_S23_L005 \
RGPL=illumina \
RGPU=WCRO84_S23_L005 \
RGSM=WCRO84_S23_L005

```

Or is a better way to determine read groups?

Thank you in advance,

# Answer
> 4 votes
Usually it doesn't actually matter, but if you want it to be very correct:

```
RGID=SomeLibraryID \
RGLB=SomeLibraryID \
RGPL=illumina \
RGPU=Some_Sequencer_ID \
RGSM=SomeSampleID

```

That assumes you have one library per sample and aren't splitting the output BAM files by lane.

---
Tags: sequence-alignment, sam, bwa, picard
---


# Question
Title: Which gene I should select from this qqplot
I have a qqplot of my whole genome sequencing data; A plot is for showing possibly significant driver genes. I tried to read about qqplot though but people only say about the skewedness while I want to know from these two genes which one are more likely to be driver

# Answer
Given that TP53 is the most significant and is already known to have driver mutations in cancer it would seem to be the logical choice.
> 3 votes

# Answer
You can't know which is the driver and which is the carrier. At most you can say that a specific gene deviate more of the expected underlying hypothesis. See also other resources online.

You also seem to ignore other genes that deviate more of your null hypothesis. I recommend to plot the histogram of the p-values to see if the distribution is uniform
> 2 votes

# Answer
This looks to be a pretty good Q-Q plot. My recollection of doing Q-Q plots is that a good Q-Q requires a "more or less" linear relationship for the model to be considered okay. Again this model looks good, its 1:1, albeit it does deviate a little (for the gene of interest). TP15 falls outside the 1:1.

@Devon Ryan states TP15 has good biological credentials, which is clearly important. It maybe worth checking MAP15 because this is "more or less" within your model (for my Q-Q plots definitely within the model :-) ).

In summary, its just my opinion - that on the basis of the Q-Q plot the driver is the not the expected candidate.
> 1 votes

---
Tags: variant-calling, mutations, wgs
---


# Question
Title: Chunk alignment in a name sorted bam for parallel processing
I have a bam file with 1 billion alignment reads of which there are 700 million unique reads. I want to split the alignments into chunks for parallel-processing. Multi-alignments of the same read should be in the same chunk. I want to use pysam to solve the problem and would be very grateful for suggestions.

# Answer
python will be to slow for this job. Here's a `awk` solution. One need to sort by read name and take track over the number of reads per chunk. If the number is reached the next chunk can only created of the read name is different.

```
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
  BEGIN { part=0; line=n }
  /^@/ {header = header$0"\n"; next;}
  { if( line>=n && $1!=last_read ) {part++; line=1; printf header > part".sam"; print $0 >> part".sam" }
    else { print $0 >> part".sam"; line++; }
    last_read = $1;
  }'

```

The output will be in sam format.

You can controll the number of reads per chunk by setting the `n` value for the `awk` command.

EDIT:

A second version that will output in `bam`:

```
$ samtools sort -n -O SAM input.bam|awk -v n=1000000 -v FS="\t" '
BEGIN { part=0; line=n }
/^@/ {header = header$0"\n"; next;}
{ if( line>=n && $1!=last_read ) {print part,line; part++; line=1;}
  print line==1 ? header""$0 : $0 | "samtools view -b -o "part".bam"
  last_read = $1;
  line++;
}

```
> 2 votes

# Answer
I don't know what the timings would be like, but the python code below will output in BAM rather than SAM, so you won't earn your PIs ire for using all that disk space, and I guess your processing code if going the be the slow bit.

```
import pysam
infile = pysam.AlignmentFile("input.bam")

chunk_size = 10000000
outfile_pattern = "output_segement%d.bam"

chunk = 0
reads_in_this_chunk = 0
old_name = None
outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

for read in infile.fetch(until_eof=True):

    if old_name != read.query_name and reads_in_this_chunk > chunk_size:
        reads_in_this_chunk = 0
        chunk += 1
        outfile.close()
        outfile = pysam.AlignmentFile(outfile_pattern % chunk, "w", template = infile)

    outfile.write(read)
    old_name = read.query_name
    reads_in_this_chunk += 1

outfile.close()

```
> 2 votes

---
Tags: python, sam, pysam, multithreading
---


title: Highly heterozygous reads mapping

text:
# Question
Title: Highly heterozygous reads mapping
I have short (67bp) Hi-C reads from a highly heterozygous organism (~15% between-haplotype SNP divergence) and I have both reference haplotypes.

I wanted to try and compare different haplotyping softwares for Hi-C reads using these reads as benchmarking dataset. When mapping the reads separately to each haplotype, I get good mapping statistics. When I map the reads to a single haplotype with all the heterozygous SNPs masked (into N's), I get very poor mapping rates.

I would like to be able to map the reads when the real haplotypes are unknown (reference is a mix of haplotypes).

I use minimap2 to map the reads with the sr preset. I tried decreasing the mismatch penalty (-B) to 1 and increasing the --score-N value, but this had no effect.

As shown in the attached IGV screenshot, coverage drops to 0 when the SNP density increases. Is it feasible to map reads with such a high heterozygosity on a single (masked) reference ? Should I use another tool ?

# Answer
> 4 votes
I believe you may be able to map your reads, but I don't know how to do that with minimap2.

I recommend running gsnap, which is more SNP tolerant and provides a number of parameters which are likely to help.

For instance, I believe that most aligners will treat 'N' characters as mismatches when you align. GSNAP has parameters to account for this.

```
--query-unk-mismatch=INT       Whether to count unknown (N) characters in the query as a mismatch
                               (0=no (default), 1=yes)
--genome-unk-mismatch=INT      Whether to count unknown (N) characters in the genome as a mismatch
                               (0=no, 1=yes (default))

```

It also has a mismatch parameter similar to the one you described for minimap2.

```
-m, --max-mismatches=FLOAT     Maximum number of mismatches allowed

```

Try running with the 'genome-unk-mismatch' parameter above (with the masked reference). That may be your best bet. There are also other parameters that may help, but this should be a good start.

---
Tags: read-mapping, haplotypes
---


# Question
Title: how to obtain only one cluster's data (.csv file) in Seurat
I was using Seurat to analysis single-cell RNA Seq. And I was interested in only one cluster by using the Seurat. Does anyone know how to achieve the cluster's data(.csv file) by using Seurat or any other method that can do this function? Thank you!

# Answer
> 1 votes
You can access data within the Seurat object using `GetAssayData`, and extract a list of cell names for the cluster you're interested in using `WhichCells`:

```
library(Seurat)

cells.use <- WhichCells(object = pbmc_small, ident = 1)
expr <- GetAssayData(object = pbmc_small, assay.type = "RNA", slot = "data")[, cells.use]
expr <- as(Class = 'matrix', object = expr)
write.csv(x = expr, file = "expression_cluster1.csv", quote = FALSE)

```

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


# Question
Title: Strategies to disentangle the effect of two genes on phenotype
I am analyzing the effect of one putative tumor suppressor (and its loss) on a specific cancer phenotype. This putative tumor suppressor is known in literature to have an effect on the growth of tumors in mouse models.

However, in the same region of this gene there's a second well-known tumor suppressor - actually a very powerful one. By looking at the data from cBio portal, there's a very strong co-occurrence of loss on both tumor suppressors in a pan-cancer analysis, that is they are almost always lost together.

Therefore, I would like to disentangle the effect of the two genes, and to analyze only the effect of my putative candidate on the phenotype. The simplest approach is to use only the samples that have a single loss (only the loss of my putative candidate, but not the loss of the confounding tumor suppressor). However, only ~10% of all samples shows such behavior, and this 10% is distributed across all cancer types, decreasing the statistical power of the analysis - since I would like to test each cancer type separately to avoid mixing incompatible data.

Which kind of strategy/analysis can be used in this case?

# Answer
Maybe you could fit a mixed model to predict phenotype, using the two genes as fixed effects and cancer type as a random effect, and ask the model for the separate effects of each gene, and the interaction between them.
> 2 votes

# Answer
Or even more simply using generalized linear modeling (GLM), where the regression weights are the genes and the accuracy of the model is measured in the residual. You also Q-Q plots to assess the fit of the model to the data.

If the residual is high and Q-Q "wobbly" transformations works wonders.

---

Machine learning such as a neural network is a very cool way to do this calculation, but IMO you need Python - although people on this board do this in R. The concept is very easy you are catagorising (keep in mind I'm not an oncologist) gene vs. cancer resistance and you want neural network, linear-based regression etc.. to make the call (you will see it visually). You then use the weights to assess the relative contribution of each gene.

The 10% issue is formally assessed, its not hard to obtain the power because its a common problem. The formal way to assess this, from memory, is you assign a "null model" and the accuracy of a ML model needs to a lot higher than the null model. I forget the details, I've never been in these situations. There's a specific technical name (which I should look up).

---

The explicit GLM / GLMM (non-ML) is the easiest approach however. If you use a Mac there is a great (relatively cheap) package for GLM, but any good software package will do. I assume GLMM is available in R (certainly Stata).
> 1 votes

---
Tags: gene
---


# Question
Title: Combat for multiplatform batch correction
I have been looking ways to use data from various micro-array platform such as agilient, rosetta/merk. Some of the established method I came across is combat, but I'm not sure if it can be done for cross platform.

A similar question here where the answer is kind of can't be or shouldn't be done.

Any suggestion how to take account of batch effects across multiple platforms.

# Answer
> 3 votes
Usually with microarrays you want to make a case/control comparison, so I am going to assume that.

Data from different array platforms is generally difficult to compare: each platform is measuring potentially a different part of the expression of the gene (different exons or different regions), each platform is likely to need different normalisation approaches, and have different dynamic range.

If you have case/control data from several different platforms, then the most robust thing to do is to analyse each platform separately, and produce a list of up and down-regulated genes according to each platform.

Then, it would be reasonably safe to combine the lists of genes in a 'meta-analysis' approach. I would expect that in a lot of cases, 'up' and 'down' regulation would be correlated across platforms, but I would not necessarily expect very strong correlation in the magnitude of changes between platforms for individual genes.

---
Tags: r, microarray, batch-effects
---


# Question
Title: Help interpreting my phylogenetic tree construction of bacterial species in same genus
I constructed a phylogenetic tree using mole blast, and it created a neighbor joining tree, so I used the alignment of those sequences with mine (1-43) in mega with TN93+G+I model but I don't understand how to interpret why my species in question (29) is so much different (or maybe novel?).

Attached is my original tree (Picture 1) and my bootstrap analysis (Picture 2):

# Answer
> 3 votes
The top phylogeny is likely mid-point rooted, but I suspect the rooting is fine. The lower diagram has been rooted on the query sequence. It is taxonomic analysis using 16S which in bacteria is considered to represent species delineation at a given percentage of divergence. I thought it was 5%, which is impossible based on the phylogram in the top picture.

Bootstrap values lower than 80% mean that the clade is not supported (monophyletic group is a better term). The delineation between 10001 and 10013 against the massive ingroup of 13 16S sequences is robustly supported by bootstraps of 95% and 99%. The divergence between them is around 2%. There is no taxonomic distinction between the 13 16S sequences based on the phylogram.

One question is whether this is full length 16S? If it is there appears to be a species designation confusion in the lower clade supporting 10001 and 10013 because it is identical to to Dermococcus species and there is no way to assess which one is valid, nor whether they are really members of the same species. The very small divergence suggesting in D. abyssis might also call this into question. 10013 is in fact a robust outgroup to 10001 and the three Dermococcus spp. but the genetic distance involved is small.

However, this is not the major result. You have recovered what appears to be a large population (if that is how you sampled it) with close genetic homology to Dermococcus strain DSM20448. Given the complexity of species lineation within this genus it is impossible to say whether this represents the same species as DSM20448 but you could call it a "species complex".

Overall the population structure (if that is how it was sampled) is heavily skewed towards a bacterial species which is hitherto received minimal attention. This population skew (if that is how it was sampled), could be reversed in different ecotopes/niches/ecoregions/patients?

In summary, I wouldn't get caught up in systematics because in this genus its a bit messy. I would focus on the population structure and these results are very clear.

---
Tags: phylogenetics, phylogeny
---


# Question
Title: Where to Download Cancer Raw Reads (fastq)?
Does anyone know where I can download cancer raw reads (fastq files), tumor and Germline for non-humans?

I wanted to make a study with human data but I don't have the control access to download raw reads.

# Answer
> 2 votes
You can try to browse SRA (sequence read archive) at https://www.ncbi.nlm.nih.gov/sra. To download fastq files from SRA, you'll need their sratoolkit: https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software.

The easy command to download fastq files by their accession number is :

```
fastq-dump --split-files accession_number

```

---
Tags: rna-seq, ngs, fastq
---


title: GEOquery errors due to list input

text:
# Question
Title: GEOquery errors due to list input
I want to extract the expression matrix from the GSE object in GEOquery.

```
library(GEOquery)

gse1 <- getGEO("GSE97168",GSEMatrix=TRUE)
datExpr = exprs (gse1[[1]])
datExpr

GSM2553017 GSM2553018 GSM2553019 GSM2553020 GSM2553021 GSM2553022 GSM2553023 GSM2553024 GSM2553025
     GSM2553026

Table(GSMList(gse1)[[1]])[1:5,]

Table(GSMList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GSMList’ for signature ‘"list"’

Table(GPLList(gse1)[[1]])[1:5,]
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘GPLList’ for signature ‘"list"’

Meta(gse1)$description
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Meta’ for signature ‘"list"’

Columns(gse1)
Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘Columns’ for signature ‘"list"’

```

Why it's not working with R version 3.5.1?

# Answer
> 1 votes
Use `GSEMatrix=FALSE`, the information you want appears to be in the soft file, rather than the matrix. Yes, this is moderately annoying since it means you need two objects to describe this particular dataset.

---
Tags: single-cell, geoquery, gse, geo
---


title: calling diploid SNVs from long reads

text:
# Question
Title: calling diploid SNVs from long reads
I'd like to call diploid SNV variants from long-reads data (~80SMRTcells PacBio).

I have generated a draft reference genome for an indivudual from a heterozygous (~4%) species (Canu+Haplomerger2).

I can use this reference for some applications:

* RNAseq mapping;
* functional annotation;
* duplicate detection;
* synteny analysis

but not others:

* analysis of compound heterozygotes;
* allele-specific expression;
* variant linkage

This is mainly why I need phasing. FYI, I tried Falcon + FalconUnzip too, but the daligner phase was just taking too much time/resources and I need to find a strategy to make it work in reasonable time (working on it).

Meanwhile, for phasing, It was suggested to use HapCUT2 \- from Canu's FAQ

> (1) "Avoid collapsing the genome" in Canu to get raw diploid assembly as complete & separate as possible (2) HaploMerger2 to get high-quality reference and alternative haploid assemblies (3) HapCUT2 or other phasing tools to get the high-quality haplotype assembly based on the reference haploid assembly.

However, HapCUT2 requires the input of a diploid SNV file. Which could be done if I had short-reads data, but yet haven't found a method for doing so with in my hands just long-reads data.

Strategies I tried to get the SNVs:

* marginPhase : method paper

+ but I had a issue in running it
* Clairvoyante : method paper

+ but I currently don't have access to suitable GPU machines to build the model - and I can't use the models provided by the authors because my organism is too different from human

Any input would be greatly appreciated!

# Answer
You may want to try the longshot tool (https://github.com/pjedge/longshot) developed for calling variants in diploid genomes from long read data.
> 2 votes

# Answer
You could try aligning your reads to the draft reference genome with for example `minimap2` and calling variants with `freebayes`. It appears there is a protocol for long reads: https://github.com/ekg/freebayes#re-genotyping-known-variants-and-calling-long-haplotypes

Disclaimer: I have not run `freebayes` myself, and it may be that this is a terrible idea. I would expect that calling SNVs (Single-nucleotide variants) from PacBio reads is very challenging. However, with enough coverage, I think it could be possible.
> 2 votes

---
Tags: vcf, snp, variant-calling, long-reads, pacbio
---


# Question
Title: Error in adding data in columns of R data frame
This is my file

```
df1 = read.csv("HSC_LSC_BLAST_karyotyope.txt",header = TRUE,sep = ",",row.names = 1)

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- rep(NA, ncol(df1))
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
#metadata$Group[seq(13,16,1)] <- 'Mono'
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'

```

> ```
> >s metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy' Error in `$<-.data.frame`(`*tmp*`, Cytogenetics, value = c("Healthy",  :
> > replacement has 4 rows, data has 23
>
> ```

Im trying to add other metadata information based on clinical information apart from giving sample labeling , im certainly doing something wrong but im not sure what it is. As my idea is the another column named Cytogenetics which would have respective information status but im getting that error.

# Answer
I know this is not really an answer, but I strongly discourage you from doing what you do, it's dangerous to impute data manually because it's prone to human mistakes - that you might mess up which columns should carry with metadata.

Actually, I think I have spotted a mistake already. This is your table using the corrected version of your code (with `metadata$Cytogenetics <- NA`, see Devon's answer):

```
>metadata
        Group Cytogenetics
HSC1      HSC      Healthy
HSC2      HSC      Healthy
HSC3      HSC      Healthy
HSC4      HSC      Healthy
Blast11 Blast           NK
Blast12 Blast           NK
Blast1  Blast           NK
Blast2  Blast           NK
Blast3  Blast           NK
Blast4  Blast           NK
Blast6  Blast           NK
Blast7  Blast           NK
Blast8  Blast           NK
Blast9  Blast           NK
LSC1    Blast           NK
LSC2      LSC           NK
LSC3      LSC           NK
LSC4      LSC           NK
LSC6      LSC           NK
LSC7      LSC           NK
LSC8      LSC           NK
Blast5    LSC     Abnormal
LSC5      LSC     Abnormal

```

I sort of think that the column `LSC1` should be in the `LSC` group and visa versa the `Blast5` should be in `Blast` group. In this case it seems that it's derived from name, you can actually write a code that will extract the group from the name.

```
df1 = read.csv("HSC_LSC_BLAST_karyotyope.txt", header = TRUE, sep = ",", row.names = 1)

metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- sub("[^[:alpha:]]+", "", (colnames(df1)))

```

However, I don't know how do you know which sample is Healthy, NK or Abnormal. I suppose you have a table where this is written, if it's the case, use that table, don't impute the data manually.
> 1 votes

# Answer
Add the `Cytogenetics` columns and then assign things to a subset of it.

```
metadata <- data.frame(row.names = colnames(df1))
metadata$Group <- NA     # you don't need to write a vector; just NA does the job
metadata$Group[seq(1,4,1)] <- 'HSC'
metadata$Group[seq(5,15,1)] <- 'Blast'
metadata$Group[seq(16,23,1)] <- 'LSC'
metadata$Cytogenetics <- NA    # initiate the column with NA values
metadata$Cytogenetics[seq(1,4,1)] <- 'Healthy'
metadata$Cytogenetics[seq(5,21,1)] <- 'NK'
metadata$Cytogenetics[seq(22,23,1)] <- 'Abnormal'

```
> 2 votes

---
Tags: r
---


# Question
Title: How do I perform a pathway related grouping of genes?
I have a list of gene names from the analysis of my mouse single seq data. I want to perform a simple grouping of the genes based on the pathways to which they belong or pathways through which they relate. Results from my online searches strongly recommend the ingenuity pathway analysis (IPA) tool. My laboratory has a license for the tool, but I need help with the workflow/steps. I have no prior experience with pathway analysis, and I have a background in R and Python programming. Can anyone also suggest any link which can help me proceed with the required IPA analysis? I will also be very grateful for any other tools (databases, R or Python packages) that can be helpful here.

# Answer
The DAVID tool is always interesting, and easy to try: https://david.ncifcrf.gov/gene2gene.jsp

You submit a list of genes, it uses several different annotation databases to annotate the gene list, and then it clusters the results and tests for statistical enrichment. You can potentially learn about functional enrichment in your gene list from several different independent annotations.
> 2 votes

---
Tags: scrnaseq, bioconductor, proteins, pathway, gene-ontology
---


# Question
Title: Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3
I downloaded the annotation of the *C. elegans* genome in GFF3 format from Ensembl. I typed the following command, hoping to get the header of the file (lines starting with `#`).

```
grep '^#' Caenorhabditis_elegans.WBcel235.95.gff3

```

This returned 46798 lines, many many more than the dozen or so header lines. After searching the file for the `#` character, I noticed many lines that contain only three hashtag (`#`) symbols.

Why are these lines there?

# Answer
There are some good answers so far, but I don't think any of them fully communicate the significance of the `###` directive. The GFF3 specification states:

> *This directive (three `#` signs in a row) indicates that all forward references to feature IDs that have been seen to this point have been resolved. After seeing this directive, a program that is processing the file serially can close off any open objects that it has created and return them, thereby allowing iterative access to the file. Otherwise, software cannot know that a feature has been fully populated by its subfeatures until the end of the file has been reached. It is recommended that complex features, such as the canonical gene, be terminated with the `###` notation.*

In other words, without the `###` directive you have to load the entire GFF3 file into memory to make sure you've inferred all feature relationships completely (as encoded by `ID` and `Parent` attributes). With the `###` directive, it's possible to process the annotations in smaller chunks, with a MUCH smaller peak memory usage.
> 4 votes

# Answer
The three `#` are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the `#`.

To `grep` just the header, which has 1 or 2 `#` from the line start on, you can use extended regular expression:

```
$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

```

Which means: "grep all lines that started with 1 or 2 `#` and are not followed by `#`"

If you like to create a new `gff3` to get rid of this separator lines you can do this:

```
$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3

```

The `-v` prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")
> 4 votes

# Answer
As others have stated, those are just there to separate the entries for easier parsing. They enable you to do nifty tricks like:

```
$ awk -v RF='###' '/Y74C9A.3/' Caenorhabditis_elegans.WBcel235.95.gff3
I   WormBase    mRNA    4116    10230   .   -   .   ID=transcript:Y74C9A.3;Parent=gene:WBGene00022277;Name=Y74C9A.3;biotype=protein_coding;transcript_id=Y74C9A.3
I   WormBase    three_prime_UTR 4116    4220    .   -   .   Parent=transcript:Y74C9A.3
I   WormBase    exon    4116    4358    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e5;constitutive=1;ensembl_end_phase=-1;ensembl_phase=0;exon_id=Y74C9A.3.e5;rank=5;version=1
I   WormBase    CDS 4221    4358    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    5195    5296    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e4;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e4;rank=4;version=1
I   WormBase    CDS 5195    5296    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    6037    6327    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e3;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e3;rank=3;version=1
I   WormBase    CDS 6037    6327    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    9727    9846    .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e2;constitutive=1;ensembl_end_phase=0;ensembl_phase=0;exon_id=Y74C9A.3.e2;rank=2;version=1
I   WormBase    CDS 9727    9846    .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    CDS 10095   10148   .   -   0   ID=CDS:Y74C9A.3;Parent=transcript:Y74C9A.3;protein_id=Y74C9A.3
I   WormBase    exon    10095   10230   .   -   .   Parent=transcript:Y74C9A.3;Name=Y74C9A.3.e1;constitutive=1;ensembl_end_phase=0;ensembl_phase=-1;exon_id=Y74C9A.3.e1;rank=1;version=1
I   WormBase    five_prime_UTR  10149   10230   .   -   .   Parent=transcript:Y74C9A.3

```

To extract the headers, you just need the lines that start with exactly 2 `#`:

```
$ grep '^##[^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942

```

Or, to also include the comments, lines that start with a `#`, then either a `#` or a `!` and then a non-`#`:

```
$ grep '^#[#!][^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06

```

Or just a `#`, any character and then a non`-#`:

```
$ grep '^#.[^#]' Caenorhabditis_elegans.WBcel235.95.gff3
##gff-version 3
##sequence-region   I 1 15072434
##sequence-region   II 1 15279421
##sequence-region   III 1 13783801
##sequence-region   IV 1 17493829
##sequence-region   MtDNA 1 13794
##sequence-region   V 1 20924180
##sequence-region   X 1 17718942
#!genome-build WormBase WBcel235
#!genome-version WBcel235
#!genome-date 2012-12
#!genome-build-accession NCBI:GCA_000002985.3
#!genebuild-last-updated 2017-06

```
> 3 votes

# Answer
Hmmm, I dunno. If I was doing this I'd want an easy way to parse the file and ### would be a good marker. Anyway just pop a space character after the #. Using sed ...

```
sed -n '/^# /p' < ./Caenorhabditis_elegans.WBcel235.95.gff3

```

That should remove the ### lines

Also `head -20` and `sed -n '/1,20/p' ./Caenorhabditis_elegans.WBcel235.95.gff3` would get around the problem with a bit of tweaking.

---

The grep verion is (removing unnecessary jargon)

```
grep '^# ' Caenorhabditis_elegans.WBcel235.95.gff3

```

This is the same command,

```
grep '^#\s' Caenorhabditis_elegans.WBcel235.95.gff3

```

The point behind the code is there is no negation operator, in sed (I assume grep is the same). In Perl you would simple \[^#\]{3} (need checking) in your regex (remove all ###). In sed (or grep) you can exclude them using the ^ (beginning of a line) and suffixed with a "space" with only a single character in the middle all 3 letter patterns are removed.

---

Another approach without the space character is

```
sed -n '/^###/d; /^#/p' ./Caenorhabditis_elegans.WBcel235.95.gff3

```

Albeit this may require,

```
sed -i '/^###/d; /^#/p' ./file

```
> -1 votes

---
Tags: ensembl, sequence-annotation, gff3, gff
---


# Question
Title: Separation of mixed plasmid DNA sequences post whole-plasmid sequencing
Imagine a DNA sample containing a mixture of different intact plasmids. These samples are sequenced using either MiSeq or HiSeq sequencing. Would it possible to assemble these plasmids post-sequencing as they would have been when sequenced individually?

# Answer
You would expect to have high coverage, given the plasmids are short, so de novo assembly would be likely very easy. Given that each plasmid is present in different multiples, you would expect different coverage on each plasmid, so it might be best to approach it as a metagenome-type or transcriptome assembly, rather than a classic genome assembly.

Alternatively, Google search just told me there is a tool called plasmidSPADES designed for just your purpose. SPADES itself is a reliable genome assembler, so I would imagine plasmidSPADES to be a good choice. It does assume that there is also some whole bacterial genomes in the sequencing mix though.

As @terdon commented, it does depend on whether you have many identical copies of each plasmid, or many similar but not identical plasmids in the mix.
> 4 votes

# Answer
Let's take a step back and consider the "perfect" output for a *de novo* assembly algorithm. Ideally, you would like to see one complete sequence for molecule (chromosome, plasmid, etc.). In reality, this is difficult to achieve due to a couple factors.

* By random chance some regions may have low coverage, meaning that there are an insufficient number of reads spanning the region to reconstruct it *de novo*.
* A bigger problem is that some portions of the genome are repetitive and occur at many different locations, sometimes even on different chromosomes.

So in practice, it's very common that *de novo* assemblies can only recover partial fragments of each molecule, and sometimes fragments are erroneously fused due to shared repetitive content.

**Back to your question**: it's certainly possible to assemble distinct plasmids that were sequenced together without any special sample prep. How well this will work in practice will depend on the amount of coverage/redundancy in the sequenced reads and the amount of DNA shared between the plasmids. If they don't share much in common, and you sample to a sufficient depth of coverage, you shouldn't have any problem assembling them.
> 2 votes

---
Tags: ngs, assembly, genome-sequencing
---