title: How to convert Bed file to fasta file?

text:
# Question
Title: How to convert Bed file to fasta file?
I have Bed file containing start and end of a sequence, and I need to convert it to fasta format, any recommendations?

# Answer
`bedtools getfasta` can do this. See documentation.
> 5 votes

# Answer
You could start with a folder of `samtools faidx`-indexed FASTA files for your genome of interest:

```
$ mkdir /genomes/hg38 && cd /genomes/hg38
$ wget ftp://hgdownload.cse.ucsc.edu/goldenPath/hg38/chromosomes/*.fa.gz
$ for fn in `ls *.fa.gz`; do gunzip $fn; done
$ for fn in `ls *.fa`; do samtools faidx $fn; done

```

Then use that folder with the following script:

```
#!/usr/bin/env perl

use strict;
use warnings;
use Getopt::Long;

#
# bed2faidxsta.pl
# --
# Reads BED data from standard input, writes FASTA to standard output.
#
# Dependent on samtools and indexed FASTA data files located in $fastaDir
# variable. Set --fastaDir=dir to set custom directory containing a source
# of per-build, bgzip-compressed FASTA and associated index (fa.gz.fai)
# files, or leave unset to use data in current working directory. Use the
# --fastaIsUncompressed option if the FASTA files are not compressed.
#

# test if samtools is available
`samtools --version` || die "Error: The samtools application is required to run this script. Try 'module add samtools' or install a local copy of samtools.\n";

# default FASTA input is current working directory
my $fastaDir = `pwd`; chomp $fastaDir;
# default is to assume input coordinates use zero-based index scheme
my $oneBased;
# default is to leave IDs alone
my $useIDPrefixAsStrand;
# default is to assume FASTA files are bgzip-compressed
my $fastaIsUncompressed;

GetOptions ('fastaDir=s' => \$fastaDir, 'oneBased' => \$oneBased, 'useIDPrefixAsStrand' => \$useIDPrefixAsStrand, 'fastaIsUncompressed' => \$fastaIsUncompressed);

if (! -d $fastaDir) { die "Error: FASTA directory does not exist\n"; }

while (<STDIN>) {
    chomp;
    my ($chr, $start, $stop, $id, $score, $strand) = split("\t", $_);
if (!defined($chr) || !defined($start) || !defined($stop)) { die "Error: No chromosome name, start or stop position defined\n"; }
    if (!defined($id)) { $id = "."; }
if (!defined($score)) { $score = "."; }
if (!defined($strand)) { $strand = "+"; } else { $strand = substr($strand, 0, 1); }
# adjust coordinates to one-based index, if necessary
my ($queryChr, $queryStart, $queryStop) = ($chr, $start, $stop);
if (!$oneBased) {
        $queryStart++;
}
# adjust strand if required
if ($useIDPrefixAsStrand) {
        $strand = substr($id, 0, 1);
    }
    # lookup
    my $queryFn = "$fastaDir/$chr.fa.gz";
if ($fastaIsUncompressed) {
        $queryFn = "$fastaDir/$chr.fa";
}
my $queryKey = "$queryChr:$queryStart-$queryStop";
my $queryResult = `samtools faidx $queryFn $queryKey`; chomp $queryResult;
# linearize result
my @lines = split("\n", $queryResult);
    my @seqs = @lines[1..(scalar @lines - 1)];
    my $seq = join("", @seqs);
# handle reverse-stranded elements
if ($strand eq "-") {
        $seq = rc_sequence($seq);
    }
    # print to standard output
    my $header = ">".join(":",($chr, $start, $stop, $id, $score, $strand));
print STDOUT $header."\n".$seq."\n";
}

sub rc_sequence {
    my $seq = shift @_;
my $reverse_complement = reverse($seq);
$reverse_complement =~ tr/ACGTacgt/TGCAtgca/;
    return $reverse_complement;
}

```

Then:

```
$ perl bed2faidxsta.pl < in.bed > out.fa

```

This approach can be useful with stranded input, nonstandard genomes, or just to do things with a bit more knowledge of inputs.
> 1 votes

# Answer
UCSC 'kent' command line tool:

```
faToTwoBit -bed=yourCoordinates.bed genomeSequence.2bit output.fa

```

Binary for this command and others available at: UCSC kent command line tools
> 1 votes

---
Tags: file-formats
---


# Question
Title: Integrative analysis of omics studies using machine learning
I would like to use public omics datasets (ChIP-seq, RNA-seq, and ATAC-seq) from different studies to do an integrative analysis as follow:

1. Normalise samples, within each type of omics, from different public datasets.
2. Convert the normalised values into a uniform scale to make the comparison between ChIP-seq, RNA-seq and ATAC-seq possible.
3. Feed the normalised uniformed values into machine learning to infer one feature (e.g. RNA expression) from other features (e.g. TF or histone marks ChIP-seq).

How could I integrate these data-sets to predict one group of variables from another one?

# Answer
> 2 votes
The steps you describe are correct. For step 2 it is usually normalized to mean 0 and variance 1. However the "machine learning" part is important.

Having several samples being technical replicates will make the integration task easier. However, you have too few samples to make any good prediction. At most I would describe it as an exploratory analysis.

I recommend using regularized canonical correlation analysis for this kind of multi-omics analysis. The method seeks which are the features of one dataset that correlate (or covariate, depending on a parameter) with the variables of other datasets. Then you could use those features to predict the position of the other.

There are several implementations of this methods, the original in the RGCCA package, and mixOmics for other procedures with cross-validations also STATegRa has other related methods.

---
Tags: rna-seq, chip-seq, atac-seq, machine-learning, multi-omics
---


# Question
Title: different results coming from biomart online and biomaRt R library
I am playing with biomart in order to get the list of genes associated to a specific gene ontology (GO) term.

So, if I use the webservice as in this query:

I get at least 50 genes. If I download the file I can gen an exact count of them:

```
$ cat mart_export.txt | wc -l
347

```

Now, I can do the same operation by using the R library biomaRt. Specifically:

```
> ensembl <- useEnsembl(biomart="ENSEMBL_MART_ENSEMBL", dataset='hsapiens_gene_ensembl')
> gene <- getBM(attributes = c('external_gene_name'), filters = 'go', values = 'GO:0030098', mart = ensembl)
> unique(gene)
  external_gene_name
1               SPI1
2              PRKDC
3                TOX
4               LY6D
5               CBFB
6              IKZF1
7               RELB
8             FLT3LG

```

By doing so, these are all the genes I retrieve. The two results are totally different and I don't know why it. All these eight genes are found in the biomart file (I grepped it, to be sure).

I believe this is not an issue related to databases versions since I checked:

```
>listMarts()
               biomart               version
1 ENSEMBL_MART_ENSEMBL      Ensembl Genes 94
2   ENSEMBL_MART_MOUSE      Mouse strains 94
3     ENSEMBL_MART_SNP  Ensembl Variation 94
4 ENSEMBL_MART_FUNCGEN Ensembl Regulation 94

```

and:

```
>listDatasets(ensembl)
...
38 Human genes (GRCh38.p12)
...

```

Both are the same version used by biomart online.

What is going on? Is it possible that biomart online is considering GO relations that are not specified in the R version?

# Answer
Sorry, this was my mistake in the last question. To search down the ontology, rather than just for the specific association with a term, biomaRt needs a different filter: 'go\_parent\_term'.

Try:

```
gene <- getBM(attributes = c('external_gene_name'), filters = 'go_parent_term', values = 'GO:0030098', mart = ensembl)

```
> 4 votes

---
Tags: r, biomart, gene-ontology
---


# Question
Title: Delete all 4 lines of a fastq read from a fastq file using read ID
I have the following error when running bowtie2:

> Error: Read HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG has more read characters than quality values.

I now want to remove all 4 lines of this specific read from the fastq file.

How can I use awk or sed to do this?

# Answer
**Don't do it:** Your FASTQ file is either malformed or your FASTQ record spans more than four lines, which is allowed in FASTQ. For a detailed description of what can go wrong in FASTQ parsing see for example http://biopython.org/DIST/docs/api/Bio.SeqIO.QualityIO-module.html#FastqGeneralIterator.

If the FASTQ is malformed, then you should really ask yourself how this happened in the first place and fix the source of the problem. If the record is valid FASTQ, then I suggest parsing the read with for example `FastqGeneralIterator` and dumping the parsed result back to FASTQ in a 4-line-per-record form.
> 9 votes

# Answer
If you are 100% sure the read only has 4 lines (they can have more), you can use this `sed` command:

```
sed -i.bak '/^@HWI-D00466:116:CC62WANXX:3:1102:7363:63646 1:N:0:GCACACG/,+3d'

```

The `-i.bak` makes `sed` modify the original file and create a backup copy with the same name and the extension `.bak`. The command just means "delete the line matching the pattern and the next three lines".
> 6 votes

# Answer
When I can I like to do file processing line-by-line in the spirit of UNIX tools. You can read 4 lines from a Fastq file into 4 tab-separated values on a single line using `paste`, and then use `grep` to filter out the record in question. (Then you just have to turn the tabs back into line breaks.)

```
paste - - - - < reads.fastq \
    | grep -v 'HWI-D00466:116:CC62WANXX:3:1102:7363:63646' \
    | tr '\t' '\n' \
    > reads-fixed.fastq

```

However, it's likely that other reads in your Fastq file are corrupted as well, in which case it's probably better to write a script in Python or some other language that discards all reads whose read length doesn't match the quality string length.

But of course the most important question is: how did the data become corrupted in the first place? Where did this file come from, and what pre-processing steps were performed on it? Is getting rid of this read (or all reads with mismatched sequence/quality length) going to solve the problem for you? These aren't questions that we can answer for you, but probably deserve good answers. :-)
> 3 votes

# Answer
Here's an `awk` one-liner that works as intended, unlike the `bioawk` answer.

This removes all instances of fastq records in which the length of the seq and qual fields do not match. Again, append `| gzip > filtered.fastq.gz` to gzip the output.

```
zcat bad_file.fq.gz | awk '{ pos = NR%4; if ( pos == 1) {h=$0} else if (pos == 2) \
  {s=$0} else if (pos==3){c=$0} else if ( pos==0 ) \
  {q=$0; if ( length(q) == length(s) ){printf("%s\n%s\n%s\n%s\n", h, s, c, q)}}}'

```

Delete the `\` characters at the end of the first two lines to make it easy to paste into your terminal.

Here is a shorter version of the same thing thanks to @terdon

```
zcat file.fastq.gz | awk '{ record[++k]=$0; if(NR%4==2){slen=length($0)} else if(NR%4==0){if(length($0) == slen){ for(i in record){print record[i]}} k=0; }}'

```
> 1 votes

# Answer
To generalize this to remove all reads that do not have the same length of read and sequence quality, you can use this bioawk one-liner:

```
bioawk -cfastx '{ if (length($seq) == length($qual) ){printf("@%s %s\n%s\n+\n%s\n", $name, $comment, $seq, $qual)} }' my_fastq.gz

```

This strips anything in the "fastq comment field", `+`, but that almost never has anything there anyway.

Append `| gzip > filtered.fastq.gz` to the above command to gzip the output.
> 0 votes

---
Tags: ngs, fastq, filtering, bioawk, text-processing
---


# Question
Title: Is it possible to use spike-ins with droplet technologies?
By reading Bacher & Kendziorski (2016), I found the following statement:

> recent droplet technologies are not yet able to accommodate spike-ins

What's the reason why current droplet technologies cannot use spike-ins? By searching through the available methodologies I found this question on 10X community, but still I cannot understand what is the problem created by using spike-ins with droplet technologies.

1. Bacher, R. & Kendziorski, C. Design and computational analysis of single-cell RNA-sequencing experiments. Genome Biol. 17, 63 (2016).

# Answer
> 1 votes
I believe it is mostly a cost issue as spike-ins are "expensive" to use given the large number of cells that are/can be sequenced via droplet based methods (as compared to plate-based ones). The 10x link that you provided actually states that this is possible with 10x and Klein et al. used ERCC spike-ins in 953 droplets (Cell. 2015 May 21;161(5):1187-1201. doi: 10.1016/j.cell.2015.04.044.).

---
Tags: scrnaseq, spike-in
---


# Question
Title: How to get PSSM of protein by protein sequence with python?
I didn't find way that convert protein sequences in FASTA to PSSM by python, and how to write program in python to find PSSM for our data set ? Maybe better to question how to find PSSM from NCBI alignments ?

Thank you

# Answer
> 1 votes
We solve our problem and create gitHub repository (bio-protein).

---
Tags: python, psi-blast, pssm
---


# Question
Title: Continuing analysis based on a log-likelihood ratio? (Comparing methylation between two PoolSeq phenotypes)
I'm trying to get some pointers on where to start on the analysis of my dataset. First, some basic info: Background: We performed Nanopore sequencing of two poolseq samples called UU and DD referring to damage vs undamaged. We were able to obtain approximately 3x coverage (this is currently more of a proof of concept than a hardcore association job, so the coverage is unfortunate but worth pursuing). Between these samples, we know that one is going to have different methylation patterns in response to a phenotypic trait response, but we don't know what regions along the genome are going to be differently methylated. The goal is to compare what's methylated (as in higher or lower levels of methylation) between the two samples, knowing there ought to be a lot that cancels out and come up with ideally a short list of significantly different regions we can associate to the plastic response.

Output of initial methylation calling data: Using the raw signal level data against the reference genome and sorted, mapped BAM files the Nanopolish software was used to call methylation. This produces a dataset with the following headers: 1. Chromosome, 2. +/- Strand 3. Start vs. 4. End of the CG dinucleotide (the difference between these is 1 or two bases) 5. The read ID the dinucleotide came from 6. the log likelihood ratio calculated from the model embedded in the methylation calling process using 7. log likelihood that it's methylated and 8. log likelihood that it's unmethylated. 8. number of calling strands (always appears to be 1) 9. number of motifs (don't know about this) and finally 10. the sequence surrounding the dinucleotide, which is always short but variable in length.

I can follow the basic interpretation of the log-likelihood... if it's a positive number, there is statistical support that it's methylated. This would probably be about as much as we need, short of a p-value, if we were looking for a specific region, but we are trying to analyze the differences in de/methylation across the entire methylome. How would you suggest using the statistic values provided in our dataset to compare and visualize the significant different modifications of the dinucleotides between strands/samples?

# Answer
I assume you don't want to spend 1-2 years developing your own method based on nanopore data, so the path forward is to convert the output into a format used by standard DMR calling programs. There's no standard format for that, annoyingly, but the general things you'll need are:

1. Coordinates of a C/CpG with coverage
2. Number of reads supporting methylation at that site.
3. Number of reads not supporting methylation at that site.

Typically you would do some filtering before/after making such a table so you only have decent quality calls (that would translate into some reasonable log likelihood for your data) and possibly only sites with decent (5x, 10x, 30x, etc.) coverage. You'd then feed that for all of the replicates in all of your groups into programs like bsseq, DSS or Metilene (to not so randomly select 3 examples).

If there aren't DMR finding programs already developed specifically for nanopore data I expect someone is developing it. The technology has a higher error rate, but given that bisulfite conversion as required in Illumina data also leads to a higher error rate that's not necessarily a huge issue if the methylation calling itself is reasonably accurate and the depth is high enough. Note that you probably need 10X coverage as a bare minimum unless your effect size is very very large.
> 2 votes

---
Tags: statistics, methylation, genetics
---


# Question
Title: What is "block-compressed" file in samtools?
SAMtools returned an error message for my gzipped genome FASTA:

> Indexed block-compressed FASTA file cannot be handled

The source code for the message is here.

```
// check if the it is a valid block-compressed file
if (IOUtil.isBlockCompressed(path, true)) {
    throw new SAMException("Indexed block-compressed FASTA file cannot be handled: " + path);
}

```

What's a "block-compressed" file?

# Answer
> 7 votes
A block compression usually refers to compressing your file into a series of small blocks (with a tool like bgzip). This allows indexing in that the index can record which record lives in which block so that the whole file does not need to be decompressed and therefore can be accessed randomly. This isn't true of files compressed as a single block using gzip.

As far as I was aware samtools only works with uncompressed FASTA files.

---
Tags: sam, samtools
---


title: How to create a .bed file from .fasta?

text:
# Question
Title: How to create a .bed file from .fasta?
I have some problems in creating a .bed file for hg19, so I will be able to visualize the .bed file in IGV.

The .fasta file contains rows of this form:

```
>chr1:0-1000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:1000-2000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:2000-3000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:3000-4000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

```

I tried both with biopython:

```
from Bio import SeqIO

with open(f1) as in_f, open(f1+'.seq','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record))

```

and awk:

```
samtools faidx file.fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' file.fasta.fai  > file.bed

```

Both ways gave me the same output which IGV couldn't work with:

```
chr1:0-1000 0   1000
chr1:1000-2000  0   1000
chr1:2000-3000  0   1000
chr1:3000-4000  0   1000
chr1:4000-5000  0   1000
chr1:5000-6000  0   1000
chr1:6000-7000  0   1000
chr1:7000-8000  0   1000

```

Now I tried to see what would be a valid .bed file format from UCSC and the format looks different. For example:

```
#bin    chrom   chromStart  chromEnd
13  chr1    39845715    39846245
19  chr1    87030756    87031903
22  chr1    111148649   111149175
22  chr1    113246065   113246628

```
1. Why my methods don't output the .bed in the same format?
2. How can I make the awk/biopython methods to output a .bed file that IGV can actually load.

# Answer
> 5 votes
Here is a fun little python script I cooked up for the occasion. It will take a standard fasta file as a command-line argument and turn it into the proper bed format. Using your example fasta:

```
import pandas as pd
import sys

inFasta = sys.argv[1] # take fasta as command argument

def fastaParser(fasta):
    headers = []
    with open(fasta) as f:
        header = None
        for line in f:
            if line.startswith('>'): # identifies fasta header line
                headers.append(line[1:-1]) # append all of the line that isnt >
                header = line[1:] # in reset header
    newHeader = (header.replace(':',',') for header in headers) # format to be accepted later
    newnewHeader = (header.replace('-',',') for header in newHeader) # format to accept later
    bedHead = (header.split(',') for header in newnewHeader) # separate by comma from format above
    return bedHead

fastaParsed = fastaParser(inFasta) # run the function

headers = list(fastaParsed) # go from generator to list
bedfile = pd.DataFrame(headers) # create dataframe which will output bedfile

bedfile.to_csv('fastabedfile.bed', sep='\t', index=False, header=None) # output bedfile format

```

to run:

```
[dkennetz@nodecn202  fasta_parser]$python fastaParser.py fakeFasta.fa

[dkennetz@nodecn202  fasta_parser]$more fastabedfile.bed
chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

```

# Answer
> 4 votes
I'm assuming the FASTA header contains all the information you want. You can easily generate a proper bed file with Biopython:

```
from Bio import SeqIO

for record in SeqIO.parse('test.fa', 'fasta'):
    lst = record.id.split(':')
    chr_, (start, end) = lst[0], lst[1].split('-')
    print('\t'.join((chr_, start, end)))

```

Note that native Python code is only marginally longer:

```
with open('test.fa') as f:
    for line in f:
        if line.startswith('>'):
            lst = line.rstrip().split(':')
            chr_ = lst[0].lstrip('>')
            start, end = lst[1].split('-')
            print('\t'.join((chr_, start, end)))

```

Both code snippets output bed format:

```
chr1    0   1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

```

# Answer
> 3 votes
The UCSC format you linked to isn't a BED file, your method should never produce it. What you posted as your desired output is also not a BED file. Below is a BED file:

```
chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

```

awk and biopython were producing the appropriate thing already. You just need to zoom in on a region in IGV. Since it's not a bigBed file, IGV won't show everything at all zoom levels.

As a demonstration, below is the BED file from above in IGV:

# Answer
> 3 votes
If your headers are **always** of the form `>chrZ:a-b`, where `chrZ` is a chromosome name, `a` is a zero-indexed, half-open start position, and `b` is a zero-indexed, half-open stop position:

```
$ awk '($0 ~ /^>/){ print substr($0,2) }' in.fa | sed 's/[:|-]/\t/g' > out.bed

```

If your headers contain coordinates using a different index, then you can make adjustments to `$1` and `$2` in the `awk` statement to create zero-indexed, half-open coordinates to export to BED, depending on the initial index.

# Answer
> 2 votes
Assuming that you have a reference genome sequence file (e.g., `reference.fasta`) available. Getting this kind of file is straightforward.

1) Index the reference genome and map your reads or FASTA sequences to it (for example with bowtie2)

```
# index reference genome (should be precomputed)
bowtie2-build reference.fasta reference
# map reads
bowtie2 -x  reference -U file.fasta -S file.sam
# compress SAM to a BAM (binary) file
samtools view -bS file.sam > file.bam

```

2) Create a BED files from the BAM file using `bedtools`

```
# extract BED file
bedtools bamtobed -i file.bam > file.bed

```

Alternatively, for sliding windows you can generate these from a reference sequence provided that you know the length of each chromosome (perhaps there is a way to extract these directly from `reference.fasta`):

```
# length per chromosome
samtools view -H file.bam | grep "SQ" | cut -d":" -f2-3 | sed 's/LN://' > file.chr.txt

```

Create sliding windows (using only the lengths per chromosome)

```
# create sliding windows
bedtools makewindows -g file.chr.txt -w 1000 > windows.bed

```

---
Tags: fasta, samtools, bed, bedtools, igv
---


title: Multiple correlation with R

text:
# Question
Title: Multiple correlation with R
I have a data frame presenting different variables in column and for different species name in rows as follow. In the column a have only one independent variables(Proteome size) and all the others are independent variables. I want to use the PIC package to correlate all these dependent variables separately to my independent variables, a way to generate the R coefficient and p-value.

```
data <- read.table("t1.txt", row.names = 1, header = T,sep='')

data
path1 path2 path3 path4 Proteome_size
ahli 4.0 -3.2 0.4 -1.0 -0.2
alayoni 3.8 3.4 -1.8 2.2 -0.3
alfaroi 3.5 2.6 0.8 -2.4 -0.5
aliniger 4.0 0.1 -1.7 2.4 0.0
allisoni 4.4 2.0 -3.7 0.5 -0.1
allogus 4.0 -2.8 0.6 -1.0 0.6
altitudinalis 3.8 2.9 -6.1 2.3 1.2
alumina 3.6 0.7 1.5 -2.7 -1.3
alutaceus 3.6 1.2 -0.6 -1.7 0.3
angusticeps 3.8 4.6 -2.0 1.2 1.7

```

I have this code following code that can correlate any one of the ‘path\*’ column to proteome size:

```
data <- read.csv("data.csv", row.names = 1
tree <- read.tree("anolis.phy")
Proteome <- data[, "Proteome_size"]
awe <- data[, "path1"]
names(Proteome_size) <- names(var1) <- rownames(data)
hPic <- pic(proteome_size, tree)
aPic <- pic(gene1, tree)
picModel <- lm(hPic ~ aPic - 1)
summary(picModel)

```

My question is that, how can I implement this code to do that for all the columns simultaneously and store the statistics info as (R and p\_value) for each variable somewhere ?

Thanks

# Answer
One option is to use one of the `apply` family functions. If you enclose your code above (which generates errors by the way) in a function, you can then "apply" this function on the rows or columns of your data with something like `apply(X, MARGIN, FUN, …)`, the margin being 1 for rows or 2 for columns (as in your case).
> 1 votes

# Answer
Like haci mentioned, you can use the apply family of functions to loop over all your columns.

```
apply(data,2, function(path_name){ summary(lm(data$Proteome_size ~ path_name -1))$adj.r.squared })

```

this gives you the adjusted r square value for your linear model, and

```
apply(data,2, function(path_name){ summary(lm(data$Proteome_size ~ path_name -1))$coefficients[,4] })

```

would give you the p-values.

What happens here is apply takes your matrix and chops it up into pieces(a piece is a column here because I gave the option 2) and each pieces is fed into the function as 'path\_name' and the output for all the pieces is given back to you.

You might get a warning about getting a perfect fit because I am correlating 'Proteome\_size' against itself.

Also,

```
awe <- data[, "path1"]
awe <- data$path1

```

are the same when dealing with columns and the second can be more convenient because you can search for autocomplete by hitting tab after typing the $
> 1 votes

---
Tags: r, correlation
---


# Question
Title: ScaleData() from Seurat causes crash on RStudio Cloud
Using the sample data in the 2,700 PBMC clustering tutorial, the session crashes at the `ScaleData()` step. The progress bar remains at 0% and hangs, followed by an R crash. I've tried reducing the size for number of genes to scale at in a single computation with the argument `block.size` with no change. The Seurat S4 object is only ~70mb so I can't imagine I'm exceeding the RStudio Cloud 1gb RAM limitations.

If you'd like to check out the cloud instance I can add you to the workspace temporarily. There is a limit on members allowed in a workspace.

I've been trying to run this tutorial on another remote cluster in RStudio inside of a container (Singularity here vs. Docker in RStudio Cloud) and am hitting the same snag.

I get receive no errors from RStudio Cloud, but in the other evirnment I get this error:

`ERROR session hadabend; LOGGED FROM: rstudio::core::Error {anonymous}::rInit(const rstudio::r::session::RInitInfo&) /home/ubuntu/rstudio/src/cpp/session/SessionMain.cpp:563`.

On this issue thread someone mentions it they receive this error message when hitting their RAM limits.

The source for `ScaleData()` can be found in the Github repo.

# Answer
> 2 votes
It seems RStudio server doesn't know it's limitations. Currently, Projects are limited to 1GB of RAM. Using the R package `ulimit` which mirrors the function of `ulimit` and allows implementation of memory limits in linux I was able to set my own limitations within the R session. After setting the limit to 1GB I see the `ulimit` limitation being tripped when invoking `ScaleData()`. Looks like a case of exceeding available system resources.

---
Tags: r, seurat, rstudio
---


title: E.coli Sequencing & Analysis

text:
# Question
Title: E.coli Sequencing & Analysis
I have been given the task of assembling a 'new' Ecoli genome and analysing the genes present etc.

The Ecoli is a new strain, and has been taken and run on a Nextseq 500 in high-output mode with 150bp paired end reads. The 'raw' files that I have is the forward and reverse reads.

I have initially QC checked the 'raw' files, and subsequently run them through trim galore and checked the QC after that.

For the next step, I now need to assemble my genome. I have been told that SPades will run a 'de novo' assembly for me, and then put that assembly into Prokka for Gene annotation.

Is this the best way to assemble the genome and annotate it? Or should I use another method? I am thinking that I should use a 'mapping' technique to assemble the genome using the Ecoli O157:H7 genome as a reference, but I have no idea how to do this. I would say that I am at an intermediate level with unix, but by no means am I a bioinformatician. Some help and guidance would be greatly appreciated!

Thanks,

Alex.

# Answer
A de novo assembly with spades should work as well as a mapping with bowtie2.

While performing a de novo assembly you should also be aware of possible contamination which you may need to remove. By using a mapping this may not be necessary.

When you use a mapping, you must be aware that you may miss relevant and diverged regions because of the mapping to a reference strain. So a mappping would be preferrable if your target genome is closely related to your reference genome.

I'd like to give you the advice to use some kind of distribution system like Anaconda which helps you to keep your system clean and make it easier to reproduce your steps.

Here are some possible steps for both ways. I am using an e.coli k12 sample, sra number SRR6940089 :

# De novo assembly

## 1. Quality Check ( you already did this )

## 2. Run Spades

```
spades.py -o k12 -m 8 -t 4 --pe1-1 SRR6940089_1.fastq --pe1-2 SRR6940089_2.fastq &> k12.log

```

This command will start an assembly using 4 threads and 8 GB RAM. You can scale it to whatever your system is capable of or what you'd like. You now have a folder called `k12` which holds the assembly information. To get some information about the contigs you can use the emboss toolkit.

```
infoseq k12/contigs.fasta

fasta::contigs.fasta:NODE_859_length_78_cov_1.000000 -              NODE_859_length_78_cov_1.000000 -              N    78     43.59
fasta::contigs.fasta:NODE_860_length_78_cov_1.000000 -              NODE_860_length_78_cov_1.000000 -              N    78     44.87
fasta::contigs.fasta:NODE_861_length_78_cov_1.000000 -              NODE_861_length_78_cov_1.000000 -              N    78     43.59

```

This is the end of the contigs.fasta. So you see contigs.fasta contains 861 contigs. And you also see thate ther length and coverage is quiet low. So you need to remove this low coverage regions.

## 3. Remove low quality contigs

There are multiple ways how you could achieve this. I'd like to filter sequence smaller than 500 bp. the tricky part is to identify the coverage threshold. Contigs with a low coverage could be some artifacts or contamination.

## 4. Check for possible contamination

So to check for possible contamination you can use kraken2. But it needs some disk space. I am using here the minikraken database

```
kraken2 --db minikraken --threads 4 --output ecoli.kraken --report ecoli_kraken.report k12/contigs.fasta &> kraken.log

```

The `ecoli_kraken.report` shows an overview of how the contigs could be classified to availible species from the database.

Based on this results you can also remove some contigs.

## 5. Run Prokka

Now after you removed unwanted contigs, you can run prokka to get the annotation.

```
prokka --cpus 4 --outdir prokka k12/contigs.fasta &> prokka.log

```

And you will find the information in the prokka folder

# Mapping

The mapping process is a little bit easier. You can use bowtie2 and samtools

```
bowtie2 -x ecoli -1 SRR6940089_1.fastq -2 SRR6940089_2.fastq -p 4 | samtools view -@ 4 -O BAM -o ecoli.bam

158216 reads; of these:
  158216 (100.00%) were paired; of these:
    35696 (22.56%) aligned concordantly 0 times
    120047 (75.88%) aligned concordantly exactly 1 time
    2473 (1.56%) aligned concordantly >1 times
    ----
    35696 pairs aligned concordantly 0 times; of these:
      4458 (12.49%) aligned discordantly 1 time
    ----
    31238 pairs aligned 0 times concordantly or discordantly; of these:
      62476 mates make up the pairs; of these:
        56512 (90.45%) aligned 0 times
        5365 (8.59%) aligned exactly 1 time
        599 (0.96%) aligned >1 times
82.14% overall alignment rate

```

You now also get some information about the similarity to your reference sequence. Now you need to work a little bit more on the BAM file. You need to sort and index this.

```
samtools sort ecoli.bam -o ecoli_sort.bam
samtools index ecoli_sort.bam

```

Now you can create a consensus sequence using bcftools and seqtk

```
bcftools mpileup -f sequence.fasta ecoli_sorted.bam | bcftools call -c | vcfutils.pl vcf2fq | seqtk seq -aQ64 -q 20 -n N > ecoli_consensus.fasta

```

This gives you now a sequence where unknown bases are marked as N.

As I wrote before, using anaconda would make the installation a lot easier. So you can write the following lines into the file `bio.yaml`

```
name: bio
channels:
  - bioconda
  - defaults
dependencies:
  - aragorn=1.2.38=h470a237_2
  - barrnap=0.9=2
  - bcftools=1.9=h4da6232_0
  - bedtools=2.27.1=he941832_2
  - blast=2.6.0=boost1.60_0
  - bowtie2=2.3.4.3=py36h2d50403_0
  - emboss=6.6.0=h6debe1e_0
  - hmmer=3.2.1=hfc679d8_0
  - infernal=1.1.2=h470a237_1
  - kraken2=2.0.7_beta=pl526h2d50403_0
  - libdeflate=1.0=h470a237_0
  - libidn=7.45.0=2
  - minced=0.3.2=0
  - parallel=20160622=1
  - perl-app-cpanminus=1.7044=pl526_1
  - perl-bioperl=1.6.924=4
  - perl-carp=1.38=pl526_1
  - perl-constant=1.33=pl526_1
  - perl-data-dumper=2.161=pl526_2
  - perl-encode=2.88=pl526_1
  - perl-exporter=5.72=pl526_1
  - perl-extutils-makemaker=7.34=pl526_2
  - perl-file-path=2.15=pl526_0
  - perl-file-temp=0.2304=pl526_2
  - perl-parent=0.236=pl526_1
  - perl-threaded=5.22.0=13
  - perl-xml-namespacesupport=1.12=pl526_0
  - perl-xml-parser=2.44=pl526h3a4f0e9_6
  - perl-xml-sax=1.00=pl526_0
  - perl-xml-sax-base=1.09=pl526_0
  - perl-xml-sax-expat=0.51=pl526_2
  - perl-xml-simple=2.25=pl526_0
  - perl-yaml=1.27=pl526_0
  - prodigal=2.6.3=1
  - prokka=1.12=3
  - samtools=1.9=h8ee4bcc_1
  - spades=3.12.0=1
  - tbl2asn=25.3=0
  - boost=1.60.0=py36_0
  - bzip2=1.0.6=h14c3975_5
  - ca-certificates=2018.03.07=0
  - certifi=2018.10.15=py36_0
  - curl=7.62.0=hbc83047_0
  - expat=2.2.6=he6710b0_0
  - fontconfig=2.13.0=h9420a91_0
  - freetype=2.9.1=h8a8886c_1
  - giflib=5.1.4=h14c3975_1
  - icu=58.2=h9c2bf20_1
  - jpeg=9b=h024ee3a_2
  - libcurl=7.62.0=h20c2e04_0
  - libedit=3.1.20170329=h6b74fdf_2
  - libffi=3.2.1=hd88cf55_4
  - libgcc=7.2.0=h69d50b8_2
  - libgcc-ng=8.2.0=hdf63c60_1
  - libgd=2.2.5=hceca4fd_3
  - libpng=1.6.35=hbc83047_0
  - libssh2=1.8.0=h1ba5d50_4
  - libstdcxx-ng=8.2.0=hdf63c60_1
  - libtiff=4.0.9=he85c1e1_2
  - libuuid=1.0.3=h1bed415_2
  - libwebp=1.0.0=h222930b_1
  - libxml2=2.9.8=h26e45fe_1
  - ncurses=6.1=hf484d3e_0
  - openjdk=8.0.152=h46b5887_1
  - openssl=1.1.1=h7b6447c_0
  - perl=5.26.2=h14c3975_0
  - pip=18.1=py36_0
  - python=3.6.7=h0371630_0
  - readline=7.0=h7b6447c_5
  - setuptools=40.5.0=py36_0
  - sqlite=3.25.2=h7b6447c_0
  - tk=8.6.8=hbc83047_0
  - wheel=0.32.2=py36_0
  - xz=5.2.4=h14c3975_4
  - zlib=1.2.11=ha838bed_2

```

Now you can run following command and most of the installation part is completed. Just the database for kraken is missing.

```
conda env create -f bio.yaml

```

I hope this gives you some help how to process your data. If you need some information just let me know ;)
> 1 votes

---
Tags: ngs, assembly
---


# Question
Title: Compare two networks of the same genes between two species
I had a set of genes for honey bees. I found the homologues of those genes in the house fly and created a network of the honey bee genes and a network of the homologous fly genes using STRING.

I know that I could do a topological comparison between the networks and check if the same set of genes give the same hubs when in two different species. However, I can't find a single study where something similar has been done.

What biological inference can I draw from the topological parameters? For example, if the diameter of the fly network is smaller than the honey bee network, what inference could I draw?

UPDATE: It would be great if atleast some publications could be given where network comparisons have been done between different species.

# Answer
So, you can give a try to this Cytoscape app that will allow to perform comparison between networks.

This app is not specifically designed for comparing networks across different species but it will enable you comparing any kind of networks.

These are the papers they published, which can give you some hints.

About your specific problem...I can't help you about what to look for when you're comparing networks from different species (not an expert in this specific topic) but as far I know you should look for something related to patterns occurring across different networks (i.e. species), network motifs, distribution of network centralities, conserved structures. In general, to compare networks, I would look for global properties (e.g. diameter, clustering coefficient...) while I would use local properties (e.g. degree, closeness, and network centralities in general) to compare the role of a node that you retrieved in all your networks.

Also, since you're looking for some literature, these links may help:
> 2 votes

---
Tags: networks
---


# Question
Title: Can I get the graph generated by cellranger
I’ve run the cellranger analysis pipeline on single cell RNASeq datasets. I can import the matrix and graph-based clusters into R. Doing this I can optimise the dimension reduction and plot cells with coloured by clusters generated by cellranger. I wish to optimise the graph-based clustering as well.

Is it possible to either:

* obtain the graph generated by cellranger from the files generated. (I can only find the final clusters in the analysis output)
* reproduce generating a similar graph (used for clustering) from the output matrix (preferably with R or Python.

With either approach, I am looking for a list of edges for nearest neighbors in the gene-barcode UMI matrix to use as input for a clustering algorithm.

# Answer
**Cell Ranger**

You can't download the tSNE coordinates for cells directly from the *Analysis* tab of the fancy, polished .html document that Cell Ranger produces. If you have access to the machine on which the pipeline was run, you can grab the intermediate files. Everything is well documented by 10X. They explain defaults etc. on that page. The tSNE coordinates are available at:

```
analysis/tsne/2_components/projection.csv

```

and the contents look like:

```
$ head -5 analysis/tsne/2_components/projection.csv
Barcode,TSNE-1,TSNE-2
AAACATACAACGAA-1,-13.5494,1.4674
AAACATACTACGCA-1,-2.7325,-10.6347
AAACCGTGTCTCGC-1,12.9590,-1.6369
AAACGCACAACCAC-1,-9.3585,-6.7300

```

**Seurat**

And the second point: Cell Ranger is replicating the Seurat workflow. Seurat uses the Louvain method for community detection by default. Check out the intro tutorial and see if you can follow along. This tutorial works with 10X data, in fact the first function call in the tutorial reads in the 10X matrices. If you are already plotting in R you'll speed through this tutorial. Seurat is great to work with, you'll have tSNEs and clusters in no time. And there's so much more you can do from there!
> 2 votes

# Answer
### Files in a cellranger run

Cellranger removes the intermediary files after the clusters have been computed. This can be observed in the following path from the cellranger project directory. This contains empty directories.

```
SC_RNA_AGGREGATOR_CS/SC_RNA_ANALYZER/RUN_GRAPH_CLUSTERING

```

However, the clusters generated by cellranger can be found in the file:

```
outs/analysis/clustering/graphclust/clusters.csv

```

Clusters can also be "exported" from `outs/cloupe.cloupe` file into a CSV using the Cloupe Browser. These CSVs can be imported and used for plotting in R or Python.

The graph structure used to compute these (with the Louvian algorithm) cannot be found in the files generated by cellranger.

### Reproducing graph-based clustering

The Suerat R package performs graph-based clustering. The shared nearest neighbour (SNN) modularity optimisation uses the same parameters as the Python scripts called by cellranger (`k=30`). This generates an SNN matrix which can be used to generate a graph structure from an adjacency matrix. For a Seurat object `test` the SNN sparse matrix can be extract with:

```
FindClusters(test, dims.use = 1:20, algorithm = 1,
             resolution = 0.6, k.param = 30, save.SNN = TRUE,
             reduction.use = "pca", metric = "euclidean")
SNN_Matrix <- as.matrix(ceiling(as.matrix(test@snn)
graph_structure <- igraph::graph_from_adjacency(SNN_matrix)

```

`FindClusters` from the `Seurat` package requires a Seurat object that has already has a dimension reduction run on it such as `RunPCA`. The `algorithm = 1` is the Louvian algorithm as performed by cellranger. The `algorithm`, `reduction.use` and `metric` can be changed as mentioned in the `Suerat` package documentation. Clustering algorithms not yet supported by Seurat can be run on the `graph_structure` resulting from extracting the SNN matrix. It is not necessary to compute the edges separately, the Seurat functions will compute and return these.

For example, the Leiden algorithm can be computed as follows:

```
devtools:install_github("TomKellyGenetics/leiden")
library(leiden)
clusters <- leiden(SNN_Matrix)

```
> 0 votes

---
Tags: scrnaseq, single-cell, clustering, 10x-genomics, cellranger
---


title: Error compiling HMMER v3.2.1

text:
# Question
Title: Error compiling HMMER v3.2.1
I tried to compile HMMER v3.2.1 in Cygwin using

```
./configure

```

followed by

```
make

```

or

```
make install

```

but failed. I got this config.log file to troubleshoot and I suspect the solution is in this line:

```
configure:2908: error: no acceptable C compiler found in $PATH

```

What does it mean? Is it really all I need to address? I followed the installation steps in this document.

# Answer
> 3 votes
I would strongly recommend to use Bioconda instead of compiling. It is one of the channels ("distributors of packages") in Anaconda Cloud which is basically a package source for conda.

In short - conda allows you to create virtual environments isolated from system and to install various precompiled software packages with one command (so, by the way, you can easily use multiple software versions on one system). Also, Bioconda is just one channel aiming on various bioinformatics software, but many more packages can be found in aforementioned Anaconda Cloud.

Tutorial how to start using Bioconda is on its homepage, and it is really easy (far more easier than to struggle with compiling).

Bioconda HMMER package

P.S. Of course there are cases in which you need a specific build of given software. Than it is really needed to compile manually.

---
Tags: software-installation, hmmer
---


# Question
Title: Is it possible to use SNP heterozygosity as a proxy for Indel heterozygosity?
I have estimated genome-wide heterozygosity levels using maximum likelihood and classical substitution model implemented in package atlas. These estimates are way more robust than classical SNP calling and simply counting number of heterogeneous locations because they integrate information from a big window into one estimate. I.e. I am quite sure about heterozygosity levels of SNPs of my genomes.

However, now I would like to know also about heterozygosity of indels and call individual snps using my robust estimates as a prior. GATK allows to specify priors on heterozygosity both of SNPs and indels, which might improve variant calls a lot, because default GATK is parametrised for human (Default priors for heterozygous SNPs is 0.001 and for indel 1.25E-4.). But, I work with insects I know that my SNPs heterozygosity is somewhere between 0.01 and 0.03 (1% - 3%).

So here comes the question, if I know the heterozygosity of SNPs, can I simply guess a reasonable prior for the heterozygosity of indels? GATK defaults for priors are 10x different, should I assume the same relation? Is there a reasonable way how to guess a prior for heterozygosity of indels?

# Answer
SNPs, insertions, and deletions are generated by different biological processes (see for example here). I would expect the per-generation rates of these mutational processes to vary across species. I would also expect the degree to which each mutational process varies across species to vary as well. I think the answer to your question is "probably not without further information".

However, priors are a feature of Bayesian statistics. In Bayesian inference, it is good practice to redo one's analysis using several priors to establish how sensitive one's results are to the chosen prior. I would suggest that you try to call indels using a range of priors and compare the quality of your call sets to determine how important the choice of prior actually is.
> 1 votes

---
Tags: snp, genomics, indel
---


# Question
Title: Using Beagle 4.1 for phasing and IBD - chromosomes
I am trying to use beagle 4.1 for phasing and IBD analyses. I was wondering if I need to run the analysis by chromosome or I can put them together in the same vcf file.

# Answer
> 4 votes
OK, I ended up contacting one of the authors of the program, who kindly answered to my question:

> Either approach should work,

\[i.e. running the analysis by chromosome vs joining the data together\]

> but I suggest dividing by chromosome because that approach has had more testing.

I hope this might help other users, I saw that several questions regarding Beagle went unanswered on biostars and seqanswers. I guess it is a pretty specialized topic.

Best,

Max

# Answer
> 1 votes
I would also add that splitting your data into chromosomes allows you to run each chromosome in parallel (presuming you have some kind of server you can send each run to), vastly reducing the amount of time your analysis is going to take. I've not run all 22 chromosomes at once, but I can imagine it will take an enormous amount of time if you are using even a small number of markers and references.

---
Tags: snp, phasing
---


title: Using Seurat to compare mutant vs.wt

text:
# Question
Title: Using Seurat to compare mutant vs.wt
I am interested in using Seurat to compare wild type vs Mutant. I don't know how to use the package. How can I test whether mutant mice, that have deleted gene, cluster together?

# Answer
> 7 votes
Single-cell analysis to compare samples is a long a difficult process. There is very good documentation for 10x Genomics cellranger, the DropSeq Pipeline and the Seurat R package. These tools all have GitHub repositories and the authors are very responsive if you encounter issues.

Depending on the technology used to generate the data, you'll need to use either cellranger or DropSeq to process the FASTQ files. These are designed to account for the the different experimental design. Still, it is important to ensure that your samples have been processed in the same way to reduce batch effects.

This question is to open-ended to cover all of the details needed to perform the analysis but here is some of the main things to consider and the general process to perform such as analysis. You will need competence using the command-line (shell) and programming in R (or Python) to perform single-cell analysis. I recommend performing this analysis on a remote server as some steps are memory intensive.

The general overview to compare samples is:

1. Demultiplex the reads for different indicies. You need to do this even if you've only put one sample per lane (to remove the index from the reads). `cellranger mkfastq` or Illumina's `bcl2fastq` will do this.
2. Obtain a reference genome (FASTA) and gene annotation (GTF) for the species you are working with. You can prepare a reference transcriptome with `cellranger mkgtf` and `cellranger mkref`.
3. Run `cellranger count` or the DropSeq pipeline on each sample separately. These will both perform STAR (splice-aware) alignment of paired-end RNA-Seq reads and count UMIs for each cell-barcode and gene. Cells will be identified and filtered by an automatic UMI threshold (this can be changed by forcing with expected number of cells). This will return a gene-barcode matrix.
4. You need to aggregate results of cellranger runs for different samples with `cellranger aggr`. This will perform downsampling by default to normalise the number of reads between samples. It will also perform graph-based (Louvian) clustering on the combined samples and return HTML and Cloupe summary data to explore, in addition to a combined gene-barcode matrix for all samples.
5. Single-cell RNA-Seq experiments are subject to batch effects. The CCA method from `Seurat` and the MNN method from `scran` are both available in the respective R packages to account for this (with different approaches). Batch effect correction needs to be performed before downstream analysis to ensure comparisons are valid.

This methods require an overlap between the samples to estimate technical errors and account for them. You need to bear in mind that neither of these methods works unless:

* some of the cells detected in each sample are expected to be the same cell type
* biological replicates of each sample have been performed

Once you have performed this correction, then standard `Seurat` workflows such as those shown in the tutorials can be performed which will identify clusters. You can use these to identify clusters specific to each genotype. Bear in mind that some overlap should be expected and if there isn't one, batch effect correction will over-correct for this.

# Answer
> 2 votes
There is no single ended answer to your question as everything will depend on your data. For example: Is this unsupervised data where you have to browse to 100 cell types for each sample or did you FACS sort on specific markers and are you looking for specific differences between 2 specific cell types etc...

Below you can find 4 recent packages that assist in your question. I strongly suggest you read the preprint of each, play around with your data and then follow the suggestion in the comments and ask more specific questions about how to proceed correctly.

If you want to use Seurat you can start here: **Seurat 3.0 pre-release to integrate your data sets:** https://satijalab.org/seurat/pancreas\_integration\_label\_transfer.html

**Seurat 2 tutorial on CCA to compare different conditions** you can use these methods on your integrated data sets. https://satijalab.org/seurat/immune\_alignment.html

**LIGER** Package and preprint can be found here https://macoskolab.github.io/liger/

**HARMONY** https://www.biorxiv.org/content/early/2018/11/05/461954

**Conos** https://github.com/hms-dbmi/conos

---
Tags: r, scrnaseq, seurat, single-cell, cellranger
---


# Question
Title: What is the minimum number of individuals per group required for ADMIXTURE analysis?
What is the minimum number of individuals required per group (population/species) for an ADMIXTURE analysis to be appropriate?

Is there a way to determine this number mathematically based on the amount of genetic data I have or other criteria?

Are there any citations I can use to figure this information out?

# Answer
The following paper provides an in-depth exploration of factors affecting power of admixture analyses under a number of different models. Although the authors state they "provide the open-source R software AdmixPower, to perform power and sample size analysis", in fact it seems the software is not available.

AdmixPower: Statistical Power and Sample Size Estimation for Mapping Genetic Loci in Admixed Populations

Tool availability:

https://research.cchmc.org/mershalab/Tools.html

Site asks for a password, apparently there is no option to sign up.
> 1 votes

---
Tags: phylogenetics
---


# Question
Title: Convert local alignments to spliced alignments in SAM file
I mapped RNA reads to reference genome, using LAST in split mode, and converted the MAF alignment to SAM with maf-convert.

My problem is that the transcripts are not reported in a spliced manner, meaning that a `transcript_ID` is reported several times in the alignment SAM file with an identical bitwise flag in `$2`. From what I understand, this is due to the fact that only the exons are mapped (one exon per row) and reported as local alignments, as the software cannot deal with splicing patterns combination exon-intron (yet), which behaviour is clear from the CIGAR string.

For a more concrete and visual example let's consider the mapping of transcript `FBtr0344900` to the reference genome as given by LAST:

```
$ cat last_aln.sam | grep FBtr0344900
FBtr0344900 0   4   42774   100 384=144H    *   0   0   TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATttttattttgatATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCAT    *   NM:i:0  AS:i:2304
FBtr0344900 0   4   43231   100 384H144=    *   0   0   CTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *   NM:i:0  AS:i:864

```

And here is the mapping of the same transcript `FBtr0344900` as given by STAR \- software which reports the alignment right the way I need it:

```
cat star_aln.sam | grep FBtr0344900
FBtr0344900    0    4    42774    255    384M73N144M    *    0    0TGCGACATTGTTCTACGATGACTACAAAAAATGACCAATAACTTCTATAAACCAATACGATATGTCAGGAGTTTCGGTCCCATACGAAGTCGCCGACTTAAGTATTTTATTTTTATTTTGATATGTGTTTGCTATTTTACCTTGTCGAATGCTTCCACACGCTATGAGAATACCATCGTGAGCGTAGCTTACTACTAGAATTTTGTTGAAGTTATTGACAAGCGATGTCTCAATATCTTCCGGACAGCCTCCAGCGTGACATTGCGGGGAATCATGTAACGGCCCAGTAACAGCCTCGGCCAGCACTCGAAGGTTTTCGTTAAGTTTAAGTATTTTATTTGTAGCACCCGCAAACAAAACATTGTGCATAAAGTCGAAGCTCATCTGGAAGCTGTTGATTGAACTGGTATTGATGGCAAGTTAAACTGGGCGACTATGTCATTTAAGGGAGATAACGCCTGAGCCGGCAGTTCTTCAATGCAGTTAACGCAATAATGCTGAGAACCGAGTATGATAATAATACACAGT    *    NH:i:1    HI:i:1    NM:i:0    MD:Z:528

```

From discussion with the author, it seems that I cannot get what I need straight from the current LAST version, and is not a technical problem. So, I have to modify the output myself. The aim would be to get at least a CIGAR line representing a complete transcript.

My question is, do you know any software to do this? What I would need is a file containing one line per uniquely mapped transcripts, featuring three fields: `transcript_ID`, `start position` and `CIGAR string`.

From my side, I proceeded like this :

1) extract the fields that I need for the interested in from the SAM file:

```
$ cut -f1,4,6
FBtr0344900  42774  384=144H
FBtr0344900  43231  384H144=

```

2) split the CIGAR line to remove items I am not interested in - I simplify the command here, just assuming I only have perfect matches (which I am interested in) and hard-clipping (which I am not interested in):

```
$ cut -f3 | sed 's/H/_ /g'| sed 's/=/= /g' | sed 's/\w*_\s*//' | sed 's/=//g'
384
144

```

3) paste the modified CIGAR in the original file, with `paste`, resulting in:

```
$ paste (1) (2) | cut -f1,2,4
FBtr0344900  42774  384
FBtr0344900  43231  144

```

4) merge rows which start ba the same `transcript_id`:

```
$ awk -F'\t' -v OFS='\t' '{ x=$1 ; $1="" ; a[x]=a[x]$0 } END { for(x in a) print x,a[x] }' |
FBtr0344900    42774    384    43231   144

```

5) compute the new cigar, calculating the intron length as the arithmetic formula `intron_length= (next_exon_start_coordinate - exon_length - previous_exon_start_coordinate)`, in this simple case above: `intron_length= 43231-384-42774`

```
$ awk '{ printf("%s", $1) }; { for (i=4; i<=NF; i+=2) { printf("%s%d%s%d", OFS, $(i-1), OFS, $i - $(i-1) - $(i-2)) } }; { printf("%s%d%s%s", OFS, $NF, OFS, RS) }'
FBtr0344900 384 73 144

```

6) ideally, with some method that I will figure out, I will modify the CIGAR string (adding alternate M,N at the end of each field except first one), this is how the final file should look like:

```
FBtr0344900 42774  384M73N144M

```

The issue of my basic approach are that:

1. I am not sure how to account for the 1-based SAM : should I add `+1` at every `exon_start_coordinate`? Does not look like, as the STAR output has the exact same cigar string that I calculated on the STAR output.
2. BIG issue: This will only work for the reads mapped in forward strand: how to make it feasible with reads mapped in the reverse strand? If I keep my current approach, I will have negative intron sizes...

Any suggestion is welcome!

# Answer
> 3 votes
Just use minimap2 in split alignment mode to realign the reads.

If that is not an option, then you could try using pysam to modify the CIGAR strings. I do not recommend this, as there are many opportunities for subtle bugs because the SAM spec is complex. You would need to:

1. Sort the BAM on read ID so that you can efficiently retrieve reads that you want to merge
2. Iterate through the BAM with pysam while loading all reads of the same read ID
3. For each set of reads with the same read ID:
    4. Order the reads of the same read ID by the number of hard-clipped bases at the start of the CIGAR string.
    5. Build a new CIGAR string from the ordered CIGAR strings
    6. Merge the meta tags from the ordered reads
    7. Write out a single, merged read

---
Tags: sam, rna-alignment, rna-splicing
---


# Question
Title: How to interpret contig-alignment.psa produced by velvet
I'm using velvet to align given reads of RNA to given CDSs (i.e. coding areas and genes) of an organism, so I can generate gene-expression profiles. But after using `velvetg out-dir/ -alignment yes`, velvet produces `contig-alignment.psa` which is a file in a strange format. The file for every contig contains zero or more records of 7 integers, and I cannot understand their meaning.

The following is an excerpt of the file, and you can see the entire file here:

```
>contig_1
1    193    1939    1939    6   194 2
195 244 1   1   6   22  71
>contig_2
1   84  1   1   6   59  142
86  170 1935    1935    6   22  106
172 285 1935    1935    6   108 221
>contig_3
1   98  2   2   6   1   98
100 321 2   2   6   100 321
334 415 1204    1204    6   1   82
>contig_15
1   23  3   3   6   84  106

```

# Answer
It looks like you are trying to do a transcript isoform expression analysis. Velvet is a genome assembler, and really not designed for this task. The tools you are looking for are kallisto and salmon.
> 2 votes

---
Tags: read-mapping, rna-alignment
---


# Question
Title: Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
After alignment using BWA, I have removed the dupliment using the samtools(Version: 1.9).

My procedure is as follows:

```
bwa mem -k 32 -M ref.fa read1 read2 > out.sam
samtools view -@ 0 -b -T ref.fa -o out.bam in.sam
samtools sort -n -o out.nameSrt.bam in.bam
samtools fixmate -r -m in.nameSrt.bam out.fixmate.bam
samtools sort -o out.fixmate.sort.bam in.fixmate.bam
samtools markdup -r -S -s in.fixmate.sort.bam out.markdup.bam
samtools flagstat in.markdup.bam > out.markdup.flagstat

```

The flagstat output result is as follows:

```
21611397 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
21611397 + 0 mapped (100.00% : N/A)
21422330 + 0 paired in sequencing
10711165 + 0 read1
10711165 + 0 read2
19797684 + 0 properly paired (92.42% : N/A)
21422330 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
1306000 + 0 with mate mapped to a different chr
727043 + 0 with mate mapped to a different chr (mapQ>=5)

```

Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output? Is there anything wrong with my procedure?

# Answer
Bwa-mem may produce chimeric alignment, where different parts of a read are mapped to distinct loci. flagstat counts them as two reads.
> 2 votes

---
Tags: ngs, samtools
---


# Question
Title: What is the most compact data structure for canonical k-mers with the fastest lookup time?
***edit: Results are current as of Dec 4, 2018 13:00 PST.***

### Background

K-mers have many uses in bioinformatics, and for this reason it would be useful to know the most RAM-efficient and fastest way to work with them programmatically. There have been questions covering what canonical k-mers are, how much RAM k-mer storage theoretically takes, but we have not yet looked at the best data structure to store and access k-mers and associated values with.

### Question

What data structure in `C++` simultaneously allows the most compact k-mer storage, a property, and the fastest lookup time? For this question I choose `C++` for speed, ease-of-implementation, and access to lower-level language features if desired. Answers in other languages are acceptable, too.

### Setup

* For benchmarking:
    + I propose to use a standard `fasta` file for everyone to use. This program, `generate-fasta.cpp`, generates two million sequences ranging in length between 29 and 300, with a peak of sequences around length 60.
    + Let's use `k=29` for the analysis, but ignore implementations that require knowledge of the k-mer size before implementation. Doing so will make the resulting data structure more amenable to downstream users who may need other sizes `k`.
    + Let's just store the most recent read that the k-mer appeared in as the property to retrieve during k-mer lookup. In most applications it is important to attach some value to each k-mer such as a taxon, its count in a dataset, et cetera.
    + If possible, use the string parser in the code below for consistency between answers.
    + The algorithm should use canonical k-mers. That is, a k-mer and its reverse complement are considered to be the same k-mer.

Here is `generate-fasta.cpp`. I used the command `g++ generate_fasta.cpp -o generate_fasta` to compile and the command `./generate_fasta > my.fasta` to run it:

```
return 0;
//generate a fasta file to count k-mers
#include <iostream>
#include <random>

char gen_base(int q){
  if (q <= 30){
    return 'A';
  } else if ((q > 30) && (q <=60) ){
    return 'T';
  } else if ((q > 60) && (q <=80) ){
    return 'C';
  } else if (q > 80){
    return 'G';
  }
  return 'N';
}

int main() {
  unsigned seed = 1;
  std::default_random_engine generator (seed);
  std::poisson_distribution<int> poisson (59);
  std::geometric_distribution<int> geo (0.05);
  std::uniform_int_distribution<int> uniform (1,100);
  int printval;
  int i=0;
  while(i<2000000){
    if (i % 2 == 0){
      printval = poisson(generator);
    } else {
      printval = geo(generator) + 29;
    }
    if (printval >= 29){
      std::cout << '>' << i << '\n';
      //std::cout << printval << '\n';
      for (int j = 0; j < printval; j++){
        std::cout << gen_base(uniform(generator));
      }
      std::cout << '\n';
      i++;
    }
  }
  return 0;
}

```

### Example

One naive implementation is to add both the observed k-mer and its reverse complement as separate k-mers. This is obviously not space efficient but should have fast lookup. This file is called `make_struct_lookup.cpp`. I used the following command to compile on my Apple laptop (OS X): `clang++ -std=c++11 -stdlib=libc++ -Wno-c++98-compat make_struct_lookup.cpp -o msl`.

```
#include <fstream>
#include <string>
#include <map>
#include <iostream>
#include <chrono>
//build the structure. measure how much RAM it consumes.
//then measure how long it takes to lookup in the data structure

#define k 29

std::string rc(std::string seq){
  std::string rc;
  for (int i = seq.length()-1; i>=0; i--){
    if (seq[i] == 'A'){
      rc.push_back('T');
    } else if (seq[i] == 'C'){
      rc.push_back('G');
    } else if (seq[i] == 'G'){
      rc.push_back('C');
    } else if (seq[i] == 'T'){
      rc.push_back('A');
    }
  }
  return rc;
}

int main(int argc, char* argv[]){
  using namespace std::chrono;
  //initialize the data structure
  std::string thisline;
  std::map<std::string, int> kmer_map;
  std::string header;
  std::string seq;
  //open the fasta file
  std::ifstream inFile;
  inFile.open(argv[1]);

  //construct the kmer-lookup structure
  int i = 0;
  high_resolution_clock::time_point t1 = high_resolution_clock::now();
  while (getline(inFile,thisline)){
    if (thisline[0] == '>'){
      header = thisline.substr(1,thisline.size());
      //std::cout << header << '\n';
    } else {
      seq = thisline;
      //now add the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        kmer_map[seq.substr(j,j+k)] = stoi(header);
        kmer_map[rc(seq.substr(j,j+k))] = stoi(header);
      }
      i++;
    }
  }
  std::cout << "  -finished " << i << " seqs.\n";
  inFile.close();
  high_resolution_clock::time_point t2 = high_resolution_clock::now();
  duration<double> time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to load the array." << '\n';

  //now lookup the kmers
  inFile.open(argv[1]);
  t1 = high_resolution_clock::now();
  int lookup;
  while (getline(inFile,thisline)){
    if (thisline[0] != '>'){
      seq = thisline;
      //now lookup the kmers
      for (int j=0; j< thisline.size() - k + 1; j++){
        lookup = kmer_map[seq.substr(j,j+k)];
      }
    }
  }
  std::cout << "  - looked at " << i << " seqs.\n";
  inFile.close();
  t2 = high_resolution_clock::now();
  time_span = duration_cast<duration<double>>(t2 - t1);
  std::cout << time_span.count() << " seconds to lookup the kmers." << '\n';

}

```

### Example output

I ran the above program with the following command to log peak RAM usage. The amount of time the lookup of all k-mers in two million sequences is reported by the program. `/usr/bin/time -l ./msl my.fasta`.

The output was:

```
 -finished 2000000 seqs.
562.864 seconds to load the array.
  - looked at 2000000 seqs.
368.734 seconds to lookup the k-mers.
     1046.94 real       942.38 user        78.96 sys
11680514048  maximum resident set size

```

So, the program used `11680514048 bytes` = `11.68GB` of RAM and it took 368.734 seconds to lookup the k-mers in two million fasta files.

### Results

Below is a plot of the results from each user's answers.

# Answer
> 6 votes
The question and the accepted answer are not about k-mer data structure at all, which I will explain in detail below. I will first answer the actual question OP intends to ask.

The simplest way to keep k-mers is to use an ordinary hash table. The performance is mostly determined by the hash table library. `std::unordered_map` in gcc/clang is one of the worst choices because for integer keys, it is very slow. Google dense, `ska::bytell_hash_map` and `ska::flat_hash_map`, `tsl::robin_map` and `absl::flat_hash_map` are much faster. There are a few libraries that focus on smaller footprint, such as google sparse and sparsepp, but those can be a few times slower.

In addition to the choice of hash table, how to construct the key is critical. For k\<=32, the right choice is to encode a k-mer with a 64-bit integer, which will be vastly better than `std::string`. Memory alignment is also important. In C/C++, as long as there is one 8-byte member in a struct, the struct will be 8-byte aligned on x86\_64 by default. Most C++ hash table libraries pack key and value in `std::pair`. If you use 64-bit keys and 32-bit values, `std::pair` will be 8-byte aligned and use 16 bytes, even though only 12 bytes are actually used – 25% of memory is wasted. In C, we can explicitly define a packed struct with `__attribute__ ((__packed__))`. In C++, probably you need to define special key types. A better way to get around memory alignment is to go down to the bit level. For read mapping, for example, we only use 15–23bp seeds. Then we have 18 (=64-23\*2) bits left unused. We can use these 18 bits to count k-mers. Such bit-level management is quite common.

The above is just basic techniques. There are a few other tricks. For example, 1) instead of using one hash table, we can use 4096 (=2^12) hash tables. Then we can store 12 bits of k-mer information into the 4096 part. This gives us invaluable 12 bits in each bucket to store extra information. This strategy also simplifies parallel k-mer insertions as with a good hash function, it is rare to insert into two tables at the same time. 2) when most k-mers are unique, the faster way to count k-mers is to put k-mers in an array and then sort it. Sorting is more cache friendly and is faster than hash table lookups. The downside is that sort counting can be memory demanding when most k-mers are highly repetitive.

---

The other answer is spending considerable (probably the majority of) time on k-mer iteration, not on hash table operations. The program loops through each position on the sequence and then each k-mer position. For an $L$-long sequence, this is an $O(kL)$ algorithm. It has worse theoretical time complexity than hash table operations, which is $O(L)$. Although hash table operations are slow due to cache misses, a factor of k=29 is quite significant. Another issue is that all programs in the question and in the other answer are compiled without `-O3`. Adding this option brings the `bytell_hash_map` lookup time from 314s to 34s on my machine.

The C program at the end of my post shows the proper way to iterate k-mers. It is an $O(L)$ algorithm with a tiny constant. The program keeps track of both forward and reverse k-mers at the same time and update them with a few bit operations at each sequence position. This echoes my previous comment "You should not reverse complement the whole k-mer". On the same machine, the program looks up k-mers in 5.5s using 792MB RAM at the peak. This 6-fold (=34/5.5) speedup mostly comes from k-mer iteration, given that the hash table library in use is known to have comparable performance to `bytell_hash_map`.

---

```
#include <stdio.h>
#include <stdint.h>
#include "khash.h"

static inline uint64_t hash_64(uint64_t key)
{ // more sophisticated hash function to reduce collisions
  key = (~key + (key << 21)); // key = (key << 21) - key - 1;
  key = key ^ key >> 24;
  key = ((key + (key << 3)) + (key << 8)); // key * 265
  key = key ^ key >> 14;
  key = ((key + (key << 2)) + (key << 4)); // key * 21
  key = key ^ key >> 28;
  key = (key + (key << 31));
  return key;
}

KHASH_INIT(64, khint64_t, int, 1, hash_64, kh_int64_hash_equal)

unsigned char seq_nt4_table[128] = { // Table to change "ACGTN" to 01234
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 0, 4, 1,  4, 4, 4, 2,  4, 4, 4, 4,  4, 4, 4, 4,
  4, 4, 4, 4,  3, 4, 4, 4,  4, 4, 4, 4,  4, 4, 4, 4
};

static uint64_t process_seq(khash_t(64) *h, int k, int len, char *seq, int is_ins)
{
  int i, l;
  uint64_t x[2], mask = (1ULL<<k*2) - 1, shift = (k - 1) * 2, tot = 0;
  for (i = l = 0, x[0] = x[1] = 0; i < len; ++i) {
    int c = (uint8_t)seq[i] < 128? seq_nt4_table[(uint8_t)seq[i]] : 4;
    if (c < 4) { // not an "N" base
      x[0] = (x[0] << 2 | c) & mask;                  // forward strand
      x[1] = x[1] >> 2 | (uint64_t)(3 - c) << shift;  // reverse strand
      if (++l >= k) { // we find a k-mer
        uint64_t y = x[0] < x[1]? x[0] : x[1];
        khint_t itr;
        if (is_ins) { // insert
          int absent;
          itr = kh_put(64, h, y, &absent);
          if (absent) kh_val(h, itr) = 0;
          tot += ++kh_val(h, itr);
        } else { // look up
          itr = kh_get(64, h, y);
          tot += itr == kh_end(h)? 0 : kh_val(h, k);
        }
      }
    } else l = 0, x[0] = x[1] = 0; // if there is an "N", restart
  }
  return tot;
}

#include <zlib.h>
#include <time.h>
#include <unistd.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
  khash_t(64) *h;
  int i, k = 29;
  while ((i = getopt(argc, argv, "k:")) >= 0)
    if (i == 'k') k = atoi(optarg);
  h = kh_init(64);
  for (i = 1; i >= 0; --i) {
    uint64_t tot = 0;
    kseq_t *ks;
    gzFile fp;
    clock_t t;
    fp = gzopen(argv[optind], "r");
    ks = kseq_init(fp);
    t = clock();
    while (kseq_read(ks) >= 0)
      tot += process_seq(h, k, ks->seq.l, ks->seq.s, i);
    fprintf(stderr, "[%d] %.3f\n", i, (double)(clock() - t) / CLOCKS_PER_SEC);
    kseq_destroy(ks);
    gzclose(fp);
  }
  kh_destroy(64, h);
  return 0;
}

```

---
Tags: k-mer, c++
---


# Question
Title: Number of reference sequences in a SAM file
From a single record, I can get the reference sequence ID from the field rID, but how can I get the total number of different reference sequences stored in the whole SAM file? It can be simple as just looping through all records, but I need some other solution.

# Answer
There are multiple ways. Here are two:

1. `samtools idxstats ‹bamfile›`
2. `samtools view -H ‹bamfile› | grep '^@SQ'`

Both commands give you a list of the references in the BAM file. To get their number, simply append `| wc -l`. But note that `idxstats` also outputs a final entry, `*`, which doesn’t correspond to any reference but to all unmapped reads. So subtract that from the number.

However, beware that while this lists all references associated with a BAM file, not every reference necessarily has records associated with it. So there might be no reads mapping to a given reference. This is the third column in the `idxstats` output. So, to get only those references with at least one mapped read, we can filter for those:

```
samtools idxstats ‹bamfile› | awk '$4 != 0' | head -n-1 | wc -l

```

Oh, and these (and other) commands work on **BAM** files, not SAM. The `idxstats` command additionally needs the BAM index, since it’s actually reading information from that index file, not the BAM (this also means that it only works on sorted BAM files). Working on SAM would require looping through the whole file, which is quite inefficient. You should generally be working with BAM rather than SAM, so this is hopefully not a restriction.

To do the same on a SAM file, you need to iterate over the file:

```
awk '{print $3}' ‹samfile› | sort -u | wc -l

```

Or:

```
awk '{refs[$3] = 1} END {for (ref in refs) print ref}' ‹samfile› | wc -l

```
> 3 votes

# Answer
A SAM file stores alignments to a set of reference sequences. However, those reference sequences are not usually stored in the file. Instead, they will be in a separate sequence file that is usually in FASTA format. You can get the number of reference sequences stored in that file using the command

```
grep -c '^>' reference.fasta

```
> 1 votes

---
Tags: reference-genome, sam
---


# Question
Title: Script to allow gene set enrichment analysis of 10x genomics data in R
I have 10x single cell RNA seq data. Which R package is best suited for analysis of the 10x data matrix. What is the script to prepare the data for downstream GSEA analysis.

I have already processed samples for single cell data (for multiple other purposes). my intention is to focus on subpopulations i.e. cell type x at different time points i.e. time 0, 1, 2 etc and compare the enriched gene pathways in this one cell type at different time points. there are three files that are used (at least in the R Seurat package) with the following file extensions: .mtx, .tsv (barcodes.tsv, genes.tsv and matrix.mtx). the R Seurat package merges these files to create a "Seurat object" which helps mainly with identifying individual diffl gene expression in each barcode cell

I found a package for R: goseq that at least has been used for bulkRNA seq but I was hoping someone has used this package specifically for single cell RNA seq data, and what script has been used to merge the 3 files to create an object for analysis. pardon my terminology, but I am a novice to R and have experience primarily through package use.

# Answer
> 5 votes
Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.

Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable. The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.

# Answer
> 2 votes
There is no purpose-built R package to perform gene set enrichment analysis on single-cell data but there does not need to be. You should be able to tools developed for bulk-RNA-Seq or microarray data, although you may not get as significant results from a sparse scRNA-Seq matrix as single-cell technologies have poor sensitivity and miss genes. What you need for a gene set enrichment analysis is:

* a database of pathways (or gene sets) to test enrichment for. Packages such as `GO.db` or `reactome.db` from Bioconductor support these in R.
* a gene expression matrix which you can get from the output of `cellranger` (as `matrix.mtx`). This should be normalised and corrected for batch effects. The `Seurat` R package is recommended to do this. Calling `object@data` will return the normalised seurat data matrix to pass it to other tools.
* an annotation to identify the genes (depending on the reference data you used for the cellranger "transcriptome").
* to define groups for differential expression analysis. You can import the clusters generated by `cellranger` in `outs/analysis` or export clusters of your choosing from the Loupe Browser as CSV files that can be imported into R. These should match the `barcodes.tsv` of columns of your data matrix to define the cluster for each cell barcode.
* to run a differential expression analysis to compare clusters. You can use clusters to define cell groups and perform analysis with tools you would use for bulk-RNA-Seq such as `limma`, `edgeR`, or `deseq`. Don't forget to adjust for multiple tests. You also need to correct for batch effects between replicates and log-transform you data to ensure the results reflect biological differences. This will give you a list of differentially expressed genes.
* to calculate gene set enrichment of differentially expressed genes in a gene set list such as those supported by pathway databases. You can do this like any other list of significant genes with tools used for bulk-RNA-Seq data, either with dedicated webtools or R packages such as `safe` or GSEA. You can also run the hypergeometric test in R yourself.

The `cellranger` output gives you all the data you need to use standard tools for which there are already guides on how to do differential analysis. These analyses do not differ from bulk-RNA-Seq analysis once the data has been normalised and the clusters have been defined.

\[edit:corrected typo in package name\]

# Answer
> -2 votes
I realized that your question is specifically gene set enrichment analysis. However, I personally have never tried any R package doing this. Although topGO may seem to be the one that may work for you? What I have been doing is that for gene enrichment analysis using Gene Ontology, simply paste and run, or DAVID when it was still updating. Sorry, I am not an expert in incorporating biological significance study in the computational result.

It seems that you are also not quite sure what tools you should use for downstream analysis. I think the following publication may be helpful.

https://www.ncbi.nlm.nih.gov/pubmed/30228881

The same primary author also published a paper titled, *Cluster Headache: Comparing Clustering Tools for 10X Single Cell Sequencing Data* (For some reason I could not find the link to it, will edit the answer if I find it. )

Hope that helps.

---
Tags: r, scrnaseq, 10x-genomics, gsea, cellranger
---


# Question
Title: Question about the dots on Quartile groups in boxplot
I have Microarray Normalized Expression data for a specific Gene. It looks like below in a dataframe `B`

```
SampleID    Gene    Type
Sample1     5.02    Tumor
Sample2     5.06    Tumor
Sample3     5.1     Tumor
Sample4     5.11    Tumor
Sample5    5.127    Normal
Sample6     5.12    Normal
Sample7    5.138    Normal
Sample8    5.149    Normal

```

I see that the minimum expression value in the table is 5.0 and maximum expression value is around 5.9.

I wanted to show the expression between two conditions with a boxplot and used following code.

```
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               add = "jitter", ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

```

This gave a plot like this

But when I remove the `jitter` from the code,

```
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

```

I see many black dots in the middle of the boxplot like this

May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?

# Answer
> 3 votes
You add the points with `geom_point()`. Just remove it and you will get your "empty" boxplot.

```
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))

```

Unfortunately I couldn't use `stat_compare_means(method = "t.test")` and `stat_n_text()`, but this two just add the labels.

---
Tags: r, rna-seq, microarray, ggplot2
---


title: What is the "primer masterfile" format?

text:
# Question
Title: What is the "primer masterfile" format?
I am working with an amplicon kit from Swift Bioscience and they provide the primer locations in what they call a "masterfile". The format is like this:

```
$ head -n2 Acel-Amplicon_56g_masterfile.txt
1       43814951        43815064        MPL_S505N-A519T_F/2     43814932        43814951        MPL_S505N-A519T_F/2_primer      43815064  43815081 MPL_S505N-A519T_R/1_primer      GCCTGGATCTCCTTGGTGA     CTTCGGCTCCACCTGGT
1       115256480       115256576       NRAS_G60-S65_new_R/1    115256457       115256480       NRAS_G60-S65_new_R/1_primer     115256576 115256602        NRAS_G60-S65_new_F/2_primer     TGATGGCAAATACACAGAGGAAG CAGGATTCTTACAGAAAACAAGTGGT

```

This is the format which is understood by Swift's primer clipping tool, but in order to use other tools, I would like to convert it to BEDPE so I can then use standard tools to clip the primers from my bam files. However, I can't do that since I am not sure what the fields in this file are.

This is what I think thy are:

1. Chromosome
2. Target region start
3. Target region end
4. 5' Primer name
5. 5' Primer start
6. 5' Primer end
7. 3' Primer name
8. 3' Primer start
9. 3' Primer end
10. 5' Primer sequence
11. 3' primer sequence

Is my interpretation correct? Also, is this a standard file format or something specific to the Swift tool? I couldn't find any documentation about a "masterfile" format online.

# Answer
> 2 votes
this is the swift standard primer file. You are a bit off in your interpretation of the columns (seems like you missed one as there are 12 columns, not 11 which threw everything off):

1) Chromosome

2) Target Start

3) Target Stop

4) **Target Name** (change which affected everything downstream in counting)

5) 5' Primer Start

6) 5' Primer Stop

7) 5' Primer Name

8) 3' Primer Start

9) 3' Primer Stop

10) 3' Primer Name

11) 5' Primer Sequence

12) 3' Primer Sequence

To convert to a BEDPE describing the primer locations, you only will need columns 1, 5, 6, 8, and 9 (and optionally 4). The last letter of the Region name for the denotes the strand orientation of the region. For example Region 1 has \_F denoting forward or (+) strand for BEDPE usage, while region 2 has \_R in name denoting reverse or (-) strand for BEDPE. Putting all this together, you can convert to BEDPE with:

```
$ awk -v OFS="\t" '{
            o=".";
            if($4~/_F\//){
                o="+"
            }
            else if($4~/_R\//){
                o="-"
            }
            print $1,$5,$6,$1,$8,$9,$4,".",o}' foo.masterfile
1   43814932    43814951    1   43815064    43815081    MPL_S505N-A519T_F/2 .   +
1   115256457   115256480   1   115256576   115256602   NRAS_G60-S65_new_R/1    .   -

```

---
Tags: file-formats, format-conversion, primer
---


# Question
Title: Double-counting coverage of overlapped read pairs
EDIT: I do not want to make any modifications to the mapped reads, I simply want to ignore one read in a read pair if they overlap the same region.

I used samtools depth to calculate the depth of coverage for samples in the whole Exome region using a GRCh37\_ref.bed. These samples are sorted and duplicate marked. I ran this calculation on a few hundred samples to determine how much more sequencing needed to be done and found something interesting. Some samples had 100X coverage at 96,000,000 mapped reads and some samples had 100X coverage at 83,000,000 mapped reads. The samples that had 100X coverage at 83,000,000 reads had read pairs overlapping certain regions of the bedfile (read 1 was covering the same coordinates as read2 to some extent).

I then found out that samtools depth double counts these overlapped regions even though they are technically from the same molecule in sequencing and would be a duplicate bases in the read. See the image for an example:

Each base in this Duplicate Overlap region is double counted when it shouldn't be. How do you guys deal with this duplicate overlap region when calculating coverage across a bam file?

samtools (v1.7):

```
samtools depth -a -b exome_targets_GRCh37.bed sample1.bam > sample1_depth.txt
awk '{s+=$3}END{print s/NR}' sample1_depth.txt

```

# Answer
`samtools depth` is a simplified version of `samtools mpileup`, which handles overlapping regions by default. Since overlapping regions shouldn't be included in coverage calculations (after all, you're not covering a region twice) you can use `samtools mpileup` instead:

```
$ cat foo.sam
@SQ SN:1    LN:1000
read1   99  1   1   50  50M =   46  95  ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0
read1   147 1   46  50  50M =   1   -95 ATTTAAAAATTAATTTAATGCTTGGCTAAATCTTAATTACATATATAATT  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<  NM:i:0

$ samtools mpileup -Q 1 foo.sam | cut -f 1,2,4
...
1   42  1
1   43  1
1   44  1
1   45  1
1   46  1
1   47  1
1   48  1
1   49  1
1   50  1
1   51  1
1   52  1
...

```

I haven't shown the entirety of the output and you'd want to use the `-aa` flag to ensure that all bases are returned.

Note that `samtools mpileup` is doing this internally by setting the base phred scores of overlapping bases in one of the mates to 0, which then get excluded due to `-Q 1` (the default is `-Q 13`, which you'd want to change).

As an aside, you probably don't need exactly correct values, only approximates. Given that, you can use `bamPEFragmentSize` from deepTools (or a similar tool from another package) to get the median fragment and read length, determine the median number of covered bases from that and use it to get a quite accurate approximation of the actual coverage (multiply the median fragment base coverage by the number of fragments and divide by the genome size). That may prove faster and close enough.
> 5 votes

# Answer
One way to deal with this would be to first merge paired-end reads based upon their overlapping regions, and then map them and calculate the coverage. This way you're only counting once per unique sequence.

Programs like SeqPrep, PEAR (Paired-End reAd mergeR), and fastq-join can do this fairly quickly. However, this isn't a perfect solution as the reads will not be joined unless they meet some minimum overlap criteria.
> 2 votes

# Answer
picard CollectHsMetrics has the parameter `PER_BASE_COVERAGE` which take this into account.

You can run it like this:

```
 $ java -jar picard.jar CollectHsMetrics \
      I=input.bam \
      O=hs_metrics.txt \
      R=reference_sequence.fasta \
      BAIT_INTERVALS=bait.interval_list \
      TARGET_INTERVALS=target.interval_list \
      PER_BASE_COVERAGE=base_coverage.csv

```

`BAIT_INTERVALS` and `TARGET_INTERVALS` could be the same file. You will need to convert your `bed` file to it using picard BedToIntervalList.

Another useful tool is mosdepth.

```
$ mosdepth --by regions.bed SampleName input.bam

```

fin swimmer
> 2 votes

# Answer
Bases from overlapping reads of the same fragment provide extra information, as their base qualities come from different sequencing cycles on the sequencer. This information can be used by downstream tools. For example, the merging tool of USEARCH increases base quality scores of bases that match between overlapping reads. Depending on the context under which coverage information will be used in the future, it may be a good idea to "double-count" overlapping bases.
> 1 votes

---
Tags: samtools, coverage
---


# Question
Title: Reconstruct a protein's 3D structure from a pair-wise distance matrix
What is a good method to reconstruct a protein's 3D structure from a pair-wise distance matrix of residues?

I think the question is self-explanatory. I know of this method but it seems that it uses contact maps, which actually use a certain cut-off. I don't have a cutoff based contact matrix but a matrix with actual distances.

# Answer
Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.

But I don't enough about the context of your problem to tell if/how it is applicable.
> -1 votes

---
Tags: pdb, protein-structure, 3d-structure
---


title: Simulate and test CNV workflow?

text:
# Question
Title: Simulate and test CNV workflow?
I'd like to evaluate a CNV project.

My aim is to evaluate if the scripts are sufficient for calling reasonable CNVs. I know they have a paper, but their scripts may be buggy... and all papers are biased...

While there are CNV simulators on Google search, none seem popular in the community. How should I approach this problem?

* Simulate CNV variants from a genome
* Simple evaluation on the called variants against my simulation?

# Answer
I cannot tell exactly from the post or the linked project what type of data these scripts are designed for, but I am going to assume that "Scripts involved in our workflow for detecting CNVs from WGS data using read depth-based methods" refers to short-read Illumina sequencing.

Simulations can be helpful, but they are usually imperfect. Nowadays we have long-read data available. I would suggest obtaining some long-read data for which short-read data is also available, doing CNV detection on both datasets independently, and comparing the results.

The problem with using long-read data is that any differences you find will be conflated with the biases of the technology. So, if you find a difference, then this approach will not tell you if the differences are due to the biases of the technologies, or due to bugs in the scripts. If you use simulations, then this conflation of biases remains.

Short of re-writing the pipeline, the only other suggestion I have is to borrow from software testing. Try constructing a few minimal examples that might also test some corner cases and see if the scripts behave exactly as expected.
> 1 votes

---
Tags: ngs, genome, cnv
---


# Question
Title: Exporting nt's differences from MEGA alignment
I have an alignment of some full genome sequences of DENV1 that we have recently sequenced from the serum of patients. I have already cleaned the sequences, built the alignments and done the SNPs calling. After that I got the consensus sequences from each sample and I used MEGA to build the alignments of the consensus sequences and to construct a tree etc.

From just looking the alignment of the consensus sequences I can see several nt s differences among the sequences.

Does anyone know if you can export from MEGA like a sort of table where you could see the nt changes for each position and sample?

# Answer
> 1 votes
You don't say which version of MEGA you are using, but this works for MEGA X on Linux:

1. Open your alignment with the "Sequence data explorer"
2. At the "Highlight" menu, select "Variable sites".
3. Then at the menu "Statistics", select "Use only highlighted sites" and select what you consider an appropriate output - I choose "Display results as CSV".
4. Now at the menu "Data" select "Export data"; at the "Sequence data export options window" if necessary select the data format again, and select "For each site" at "Writing site numbers", and select "Only highlighted sites" at "Selected sites to include".

This will result in a comma-separated values file, showing only positions with variable sites.

---
Tags: sequence-alignment, software-recommendation
---


# Question
Title: Best practices for deciding if two structural variants are actually the same variant?
I know that I can use `bcftools isec` to compare two VCF files containing single nucleotide variants (SNVs); it will generate four possible output files: one with all of the unique called SNVs in file 1, one with all of the unique called SNVs in file 2, one with all of the SNVs from file 1 that were also found in file 2, and finally one with all of the SNVs from file 2 that were also found in file 1.

Is there a similar type of tool that will allow one to compare structural variants (SVs) between two VCF files? I am about to tweak one of my Perl scripts (that currently uses the VCF.pm Perl module to parse VCF files) to perform this calculation (because thus far I have not found any tools designed for SVs) but I thought I would ask in this forum in case anyone else has already invented this wheel.

# Answer
The tool that you want to use for this task is `bedtools window`.
> 1 votes

---
Tags: vcf, ngs, software-recommendation
---