# Question
Title: cn.MOPS fails with 'missing value where TRUE/FALSE needed'
I am trying to use cn.MOPS to call CNVs on a set of whole exome sequencing data bam files. My script (up to the point where it fails) is:

```
source("https://bioconductor.org/biocLite.R")
library(cn.mops)
BAMFiles <- c('sample1.bam', 'sample2.bam', 'sample3.bam', 'sample4.bam', 'sample5.bam', 'sample6.bam', 'sample7.bam')
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)
resCNMOPS <- exomecn.mops(X)

```

The last command, `resCNMOPS <- exomecn.mops(X)` fails with:

```
Error in if (all(segMedianT == 0)) { :
  missing value where TRUE/FALSE needed
In addition: Warning message:
In normalizeChromosomes(X, chr = chr, normType = normType, qu = normQu,  :
  Normalization for reference sequence  chrUn_gl000228 not applicable, because of low number of segments

```

If I run `traceback` immediately afterwards, I get:

```
> traceback()
4: FUN(newX[, i], ...)
3: apply(sINI[chrIdx, , drop = FALSE], 2, segment, minSeg = minWidth,
       ...)
2: cn.mops(input = input, I = I, classes = classes, priorImpact = priorImpact,
       cyc = cyc, parallel = parallel, normType = normType, normQu = normQu,
       norm = norm, sizeFactor = sizeFactor, quSizeFactor = quSizeFactor,
       lowerThreshold = lowerThreshold, upperThreshold = upperThreshold,
       minWidth = minWidth, segAlgorithm = segAlgorithm, minReadCount = minReadCount,
       useMedian = useMedian, returnPosterior = returnPosterior,
       ...)
1: exomecn.mops(X)

```

I am using R version `3.4.4` on an Ubuntu 16.04.5 LTS machine.

My bed file only has one region which is on `chrUn_gl000228`:

```
$ grep chrUn_gl000228 region.bed
chrUn_gl000228  112604  113879

```

If I remove that region and rerun the same code, it works as expected. So, my question is, how can I detect this and fix it automatically? Given a collection of bam files and the bed file describing the capture regions used to sequence the samples and generate the bam files, how can I check if cn.MOPS will complain about any of the segments?

Presumably, once I have built these two objects:

```
segments <- read.table("region.bed",sep="\t",as.is=TRUE)
gr <- GRanges(segments[,1],IRanges(segments[,2],segments[,3]))
X <- getSegmentReadCountsFromBAM(BAMFiles,GR=gr,parallel=11)

```

I should be able to then check them and see if they have the right data and, if not, remove the offending segment from the `GRanges` object. But how can I do this?

# Answer
> 2 votes
I managed to track down what was producing this error. The problem was that my bed file only had one segment for chromosome `chrUn_gl000228`. However, the cn.MOPS method relies on normalizing read counts across multiple segments for each chromosome. Since this chromosome only had a single segment, the `exomecn.mops(X)` was failing.

So, either remove the single segment from the bed file, or use a bed file that has at least 2 segments on each chromosome. I created a fake segment on `chrUn_gl000228` and confirmed that the program worked as expected and didn't throw any errors.

---
Tags: r, bioconductor, cnv
---


# Question
Title: bedtools unable to open file despite being tab-delimited
I got the following error when attempting to intersect two files. The file below is supposed to be the -a file.

```
Error: Unable to open file compare_me_to.txt. Exiting.

```

When I attempted to make it tab-delimited through Ubuntu (see code below)

```
perl -p -i -e 's/ /\t/g' compare_me_to.txt

```

I got the following error:

```
Can't open compare_me_to: Permission denied.

```

Additionally, the ownership of the folder the files are in is of my username. I can't tell if this is a bedtools problem or a Ubuntu problem or something else.

# Answer
I suspect the reason is the "-i". The following will work,

```
perl -pe 's/  */\t/g' compare_me_to.txt > compare_me_to.txt2

```

-i will want to write onto the file and create a backup. I suspect you'll have spotted `chmod +w compare_me_to.txt`.

Just make sure that a double space is inserted before the * character, otherwise it will go wrong. Copying and pasting will avoid this.

---

This is a clearer solution because it avoids missing a space character,

```
perl -e 's/[ ][ ]*/\t/g' compare_me_to.txt > compare_me_to.txt2

```

and the solution you used .. (I did on OS X), what I'd recommend is,

```
perl -p -e 's/[ ][ ]*/\t/g' compare_me_to.txt # to see if it works by printing on screen without altering the file

perl -i -e 's/[ ][ ]*/\t/g' compare_me_to.txt

```

The -i flag will alter the original file and store a backup without requiring renaming the file.

**What's it doing?** The regex \[ \]\[ \]\* structure is saying look a space and if one or more spaces occur adjacent to each other replace everything with a single tab. \[ \]\* means a "space" that is present or absent and if present it can be present an infinite number of times. If the \[ \]\* is by itself the command risks saying replace '' (zip) or 'space' with a tab, so tabs can appear at every '' (zilch), so tabs can appear everywhere. The "\*" symbol in regex work differently to the "\*" wildcard character in shell.
> 1 votes

---
Tags: bedtools
---


# Question
Title: Tools to do VCF to MAF and MAF to VCF conversion?
Normally, I would use the vcf2maf scripts to convert a VCF to a MAF (or vice versa).

This is great software, but on my system, perl scripts with dependencies are easy to break. (Here it uses VEP.)

Are there any other alternatives to this?

# Answer
> 5 votes
Check this python script vcf2maf.py

# Answer
> 6 votes
You can use annovar to annotate the vcf, then convert it to maf using the function `annovarToMaf` of maftools bioconductor package.

---
Tags: vcf, variant-calling, maf
---


title: Assembly by stringtie

text:
# Question
Title: Assembly by stringtie
I run this cmd

```
./stringtie G1_sorted.bam -B -o G1.gtf -G Triticum_aestivum.IWGSC.42.gtf -p 4  -C G1.refs.gtf -A G1.abund.tab

```

Error is:

```
-WARNING: no reference transcripts were found for the genomic sequences where reads were mapped!
Please make sure the -G annotation file uses the same naming convention for the genome sequences.

```

How can I resolve it?

**Update:**

I have downloaded both (GTF and Refernce Genome)files from Ensemble. Can you confirm me that whether theu have same chromosomes naming? So that i can perform my mapping task of GTF File with Reference Genome which usually i remain fail to do due to difference in chromosomes naming. Here The Reference genome file is shown as:

```
 >1A dna:chromosome chromosome:IWGSC:1A:1:594102056:1 REF
CTAAACCCTAAACCCTAAACCCTAAACCCTAAACCCCTAACCCTAAACCCTAACCAAAAC

```

Where as The annotation file (gtf file ) is as

```
    #!genome-build IWGSC
#!genome-version IWGSC
#!genome-date 2018-07
#!genome-build-accession GCA_900519105.1
3B      IWGSC   gene    212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; gene_source "IWGSC"; gene_biotype "protein_coding";
3B      IWGSC   transcript      212892  214491  .       -       .       gene_id "TraesCS3B02G000100"; transcript_id "TraesCS3B02G000100.1"; gene_source "IWGSC"; gene_biotype "protein_coding"; transcript_source "IWGSC"; transcript_biotype "protein_coding";

```

Can it be confirmend that these both files are same ? If not same with chromosome naming convention , then please suggest me how can I change it ?

# Answer
> 3 votes
This seems to be a common problem. As @llrs pointed out, Ensembl, UCSC(genome browser and associated data) and NCBI records are not directly interchangeable. The ambiguity refers to annotation as well as the reference genome. The reference genome is identical, except for the chromosome naming convention, meaning a reference from Ensembl might not contain CHR prefix, whereas it will be present in one downloaded from UCSC.

To prevent such conflicts in downstream tools, it is good to stick to the reference and annotation from the same repository throughout the analysis and map the genes to gene names or whatever towards the end of the investigation.

You can download all the associated data for *T. aestivum* which I beleive is the organism you want to study from Ensembl plants FTP download.

---
Tags: rna-seq, assembly, transcriptome, gtf
---


# Question
Title: How can I assign a genomic region into a window using R?
I have asked this question at BioStars but have not gotten any suggestions so far. So, I am asking the community here, perhaps someone here has an idea how I can solve my problem.

I am dealing with a bunch of inversions and running some analysis on each of them. One of them is located on chromosome one, the starting position is 6199548 and ending position is 9699260. So the length of this variant is 3499712 bp. I want to generate non-overlapping windows based on the length of these focal positions (3499712) and assign the focal position to one of these windows that overlap with the starting and ending position of my focal position. I am achieving this by using the `GenomicRanges::tileGenome` function in `R`:

```
seqlength <-c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)
bin_length<-c(3499712)
 bins<-tileGenome(seqlength,tilewidth=bin_length,cut.last.tile.in.chrom
        = T)
            #bins_u<-unlist(bins)
            write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names
        = FALSE, row.names = TRUE) }

```

This generates non-overlapping 3499712 bps windows across the genome.

```
window_id chromosome start end length
"1" "chr1" 1 3499712 3499712
"2" "chr1" 3499713 6999424 3499712
"3" "chr1" 6999425 10499136 3499712
"4" "chr1" 10499137 13998848 3499712
"5" "chr1" 13998849 17498560 3499712

```

However, I have a major problem here that I have not been able to sort out. I want one of these windows to overlaps exactly with the starting and ending position of my focal position (6199548:9699260). Only the focal position not any other position outside this region. This means 2 windows (one before and one after the focal position) will have a different length which my analysis is fine with. I just do not know how to do it?

This is my desired output:

```
window_id chromosome start end length
"1" "chr1" 1 3499712 3499712 "*"
"2" "chr1" 3499713 6199547 2699834 "*"
"3" "chr1" 6199548  9699260 3499712 "*"
"4" "chr1" 9699261 13998848 4299587 "*"
"5" "chr1" 13998849 17498560 3499712 "*"

```

As you see, the the second window is now smaller, the third window accommodates the focal position and the forth window is a little larger than the bin size. Does anyone have any idea how can I do this task in R?

# Answer
> 5 votes
I'm not sure how you did your calculation but for me it looks like your region is **3499713** bp in length.

```
site <- GRanges("chr1:6199548-9699260")

seqlength <- c(chr1=159217232,chr2=184765313,chr3=182121094,chr4=221926752,chr5=187413619)

bin_length <- width(site)

## check length of region
print(bin_length)
##3499713

bins <- tileGenome(seqlength, tilewidth=bin_length, cut.last.tile.in.chrom=TRUE)

## extract overlapping bins
flank.bins <- subsetByOverlaps(bins,site)

## find new coordinates excluding region of interest
new.flank.bins <- setdiff(reduce(flank.bins),site)

## combine new bins with your region of interest
new.bins <- c(new.flank.bins,site)

## remove old bins
bins <- subsetByOverlaps(bins,new.bins,invert=T)

## add new bins
bins <- sort(c(bins,new.bins))

## check that it makes sense
##table(width(bins[seqnames(bins) == 'chr1']))

write.table(bins,file=paste("bins_good_",bin_length,sep=""),col.names= FALSE, row.names = TRUE)

```

---
Tags: r
---


# Question
Title: HISAT2 - are all multi mapped reads reported by default?
While going through the HISAT2 manual, one can see this:

> Two alignments for the same individual read are "distinct" if they map the same read to different places. Specifically, we say that two alignments are distinct if there are no alignment positions where a particular read offset is aligned opposite a particular reference offset in both alignments with the same orientation. E.g. if the first alignment is in the forward orientation and aligns the read character at read offset 10 to the reference character at chromosome 3, offset 3,445,245, and the second alignment is also in the forward orientation and also aligns the read character at read offset 10 to the reference character at chromosome 3, offset 3,445,245, they are not distinct alignments.
>
> Two alignments for the same pair are distinct if either the mate 1s in the two paired-end alignments are distinct or the mate 2s in the two alignments are distinct or both.
>
> Default mode: search for one or more alignments, report each

and then there is -k parameter (analogous to bowtie2 -K parameter)

which defines:

> HISAT2 searches for up to N distinct, primary alignments for each read, where N equals the integer specified with the -k parameter. Primary alignments mean alignments whose alignment score is equal or higher than any other alignments. It is possible that multiple distinct alignments have the same score. That is, if -k 2 is specified, HISAT2 will search for at most 2 distinct alignments. The alignment score for a paired-end alignment equals the sum of the alignment scores of the individual mates. Each reported read or pair alignment beyond the first has the SAM 'secondary' bit (which equals 256) set in its FLAGS field. See the SAM specification for details.

This is a bit confusing, as the first quote means that HISAT2 finds all distinct valid alignment and reports them by default, and the second quote says that when -k option is used, it will only report a subset of alignments which it would have reported by default. Am I correct in this inference or does one needs to use -k option to enable "multi-mapping"?

# Answer
You said:

> ...the first quote means that HISAT2 finds all distinct valid alignment and reports them by default.

Not exactly. The first section you quoted doesn't say it will find and report *all* distinct alignments, but rather that it would find and report one or more distinct alignments. It's a subtle difference, but an important one.

That section describes how distinct mappings are defined, not how many are reported. The second section you quoted, as well as the documentation for the `-k` parameter, are all consistent with this interpretation.

> ```
> -k <int>        It searches for at most <int> distinct, primary alignments for
>                 each read. Primary alignments mean alignments whose alignment
>                 score is equal or higher than any other alignments. The search
>                 terminates when it can't find more distinct valid alignments, or
>                 when it finds <int>, whichever happens first. The alignment
>                 score for a paired-end alignment equals the sum of the alignment
>                 scores of the individual mates. Each reported read or pair
>                 alignment beyond the first has the SAM 'secondary' bit (which
>                 equals 256) set in its FLAGS field. For reads that have more
>                 than <int> distinct, valid alignments, hisat2 does not guarantee
>                 that the <int> alignments reported are the best possible in
>                 terms of alignment score. Default: 5 (HFM) or 10 (HGFM)
>
>                 Note: HISAT2 is not designed with large values for -k in mind,
>                 and when aligning reads to long, repetitive genomes large -k can
>                 be very, very slow.
>
> ```

By default, 5 or 10 distinct mappings are reported (depending on the mode), but you can control this with the `-k` parameter.
> 3 votes

---
Tags: read-mapping
---


# Question
Title: Demultiplex nanopore reads with custom barcodes
We have a problem trying to demultiplex MinION sequences with custom barcodes. Do you have any software recommendations we can try for demultiplexing or how to demultiplex these custom barcodes with Albacore? We have tried using albacore but it only recognizes barcodes provided by Nanopore Technologies.

# Answer
> 5 votes
I've written some Perl code for demultiplexing nanopore reads using LAST. A workflow demonstrating the use of this script (and how it works under the hood) can be found here:

https://dx.doi.org/10.17504/protocols.io.14egnxw4zl5d

You feed the script a barcode file and input fastq file, and it creates a trained mismatch matrix, aligns to the barcodes, then demultiplexes the reads into separate files for each barcode:

```
./fastq-dental.pl -barcode barcode_full_PBK004.fa -mat bc.mat reads_all.fastq.gz

```

Example output:

```
         9 BC01
        39 BC02
         1 BC02_BC07
        16 BC03
        38 BC04
         2 BC04_BC06
         1 BC04_BC07
         1 BC04_BC09
       643 BC05
         2 BC05_BC07
    493175 BC06
       115 BC06_BC06
         2 BC06_BC07
    663161 BC07
         1 BC07_BC05
       117 BC07_BC07
         1 BC07_BC09
       106 BC08
         1 BC08_BC06
        75 BC09
        72 BC10
         1 BC10_BC07
        41 BC11
        44 BC12
         1 BC12_BC06
     13535 BCchim
    237796 BCnoadapt
       152 BCnone
        90 RB12A

```

While the approach designed around using the ONT barcodes, the barcode fasta file can be replaced by another file \[using barcode ID in the header\] without change the protocol.

FWIW, there is also now some support for detecting strand direction (e.g. using the strand switch primer provided by ONT), but I haven't yet got it working properly doing both demultiplexing and strand orientation at the same time.

# Answer
> 3 votes
There are two other ways in which you could demultiplex reads if guppy doesn't have the config file for the kit you are using.

1. Minibar: This works for dual index barcodes specifically
2. Porechop: This is intended for the ONT sets, but you can also replace them with your own. However, it is designed to work with single index barcodes. That means, you will have to run it once with your forward indices, then take that output and run again with the reverse, and then do some gymnastics to combine it all together

---
Tags: nanopore, minion, albacore, demultiplexing
---


title: Including Picard tool in galaxy

text:
# Question
Title: Including Picard tool in galaxy
I am running the Picard tool CollectSequencingArtifactMetrics in galaxy (locally in my pc).

The output files are saved into disk and not passed to *stdout* and (I guess) they can't be passed as output from picad to galaxy .

In the xml file I am including:

```
<outputs>
    <data name='preSummary' format="tabular"  from_work_dir="pref.pre_adapter_summary_metrics"/>
</outputs>

```

When I click on *View details*, in *"stderr"*, I see:

> \[Mon Mar 11 13:43:07 EDT 2019\] picard.analysis.artifacts.CollectSequencingArtifactMetrics done. Elapsed time: 0.19 minutes. Runtime.totalMemory()=3323461632

1 - It looks like it worked but shouldn't I see the highlighted information in *"stdout"* rather than in the *"stderr"*?
2- If it worked, How can I get my output files? I tracked the working directory printing `pwd` from the `<command>` but I guess all files are deleted immediately once the picard process has finished.

Any idea of what am I doing wrong?

What basic knowledge am I missing here?

Thank you!

# Answer
1. Picard, like most tools, prints logs to stderr.
2. It should be in the Galaxy history. Are you using a version of picard from the toolshed (if so, which one)? You might set `cleanup_job` to `never` in the config, which will prevent galaxy from cleaning the working directory. This will enable seeing what the job actually produced.

Most likely there's more in the log file that you're not showing is us and picard itself never ran to completion. In lieu of that information, what I wrote above is usually sufficient to debug things.
> 1 votes

---
Tags: picard, galaxy
---


# Question
Title: Does UCSC definition for telomere making sense?
UCSC provides a database of telomere regions for each chromosome:

> ftp://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/gap.txt.gz

It reported the telomere section in chr21 as:

```
grep chr21 gap.txt | grep telomere
585 chr21   0   10000   1   N   10000   telomere    no
952 chr21   48119895    48129895    515 N   10000   telomere    no

```

Telomere is the first 10000 bases in chr21?? But ... the region are all "N" in the reference genome!

How does UCSC confidently give out telomere regions if they don't even know anything about the sequences?

# Answer
It seems so, yes.

UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.

Taken from Wikipedia:

A telomere is a region of repetitive nucleotide sequences at each end of a chromosome, which protects the end of the chromosome from deterioration or from fusion with neighboring chromosomes. **For vertebrates**, the sequence of nucleotides in telomeres is AGGGTT, with the complementary DNA strand being TCCCAA, with a single-stranded TTAGGG overhang. This sequence of **TTAGGG is repeated approximately 2,500 times in humans. In humans, average telomere length declines from about 11 kilobases at birth\[ to less than 4 kilobases in old age**, with the average rate of decline being greater in men than in women.

2,500 x 6 (TCCCAA) = 15,000
(approx. 10K)

Two interesting additional facts:

The minimum telomere length needed to ensure human telomere protective stability in white blood cells is 3.81 kb
Mouse has a 100Kb telomere representation in USCS
> 3 votes

---
Tags: genomics
---


title: Running htseq-count over BAM files

text:
# Question
Title: Running htseq-count over BAM files
I am trying to derive an expression matrix from BAM files using `htseq-count`. These are bulk RNASeq BAM's by the way.

I have read the `htseq-count` documentation as well as `samtools` and figured that the following command *should* work:

```
samtools sort -on file.bam | samtools view -h - | htseq-count --type=exon --idattr=gene_id --additional-attr=gene_name - Homo_sapiens.GRCh38.77.gtf > counts.txt

```

The `.gtf` file here is from ensembl.

I have tried many variations of the above command but have failed to generate a count matrix. The *error* message states:

"0 SAM alignments processed"

Here are the first few lines of my .gtf file:

```
#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1       pseudogene      gene    11869   14412   .       +       .       gene_id "ENSG00000223972"; gene_name "DDX11L1";
gene_source "ensembl_havana"; gene_biotype "pseudogene";
1       processed_transcript    transcript      11869   14409   .       +       .       gene_id "ENSG00000223972"; trans
cript_id "ENST00000456328"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; transcript_nam
e "DDX11L1-002"; transcript_source "havana";
1       processed_transcript    exon    11869   12227   .       +       .       gene_id "ENSG00000223972"; transcript_id
 "ENST00000456328"; exon_number "1"; gene_name "DDX11L1"; gene_source "ensembl_havana"; gene_biotype "pseudogene"; trans
cript_name "DDX11L1-002"; transcript_source "havana"; exon_id "ENSE00002234944";

```

And all the genes in the text file have an expression value of 0. What exactly is going wrong in this step? Is the conversion to a SAM file the problem? I'm completely lost so any help would me much appreciated.

Thanks!

EDIT:: For anybody who sees this post cause they have a similar problem. Just note that the newer versions of `htseq-count` don't require sorted .bam files and, so following the editing of the `.gtf` file, all I needed to do was convert it to `.sam` using `samtools view -h` and then pipe this to `htseq-count`. Sorting the files prior to this conversion would cause errors.

# Answer
If your BAM file is using UCSC chromosome names (chr1, chr2, etc.) and your GTF file is using Ensembl chromosome names (1, 2, etc.) then htseq-count won't know how to match your alignments to any of the genes. The simplest solution is to not use htseq-count, but instead to use featureCounts, which will handle the chromosome name difference transparently (it's one of the few programs that does) and is much faster anyway. In the future, download your GTF and fasta files from the same source, don't mix and match them.
> 2 votes

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


# Question
Title: Searching tool to calculate phase/switch error rate
I'm looking for a tool which, given a truth vcf file and a test vcf file, calculates the phase/switch error rate. I performed phasing of a vcf using WhatsHap and want to compare the outcome to some ground truth phased vcf I have. I can't find a tool, and I don't want to write it myself. It seems like a rather common thing... so I'd expect it to exist, but my google-fu is weak today.

Example of test lines:

```
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample

chr10   100588  .       G       A       .       PASS    AF=0.18;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.08;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=11907     GT:PS   1|0:48005
chr10   102385  .       A       G       .       PASS    AF=0.07;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.1;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22665      GT:PS   0|1:48005
chr10   105170  .       G       A       .       PASS    AF=0.03;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.0;AMR_AF=0.19;SAS_AF=0.0;VT=SNP;DP=17140       GT:PS   0|1:48005
chr10   105365  .       G       T       .       PASS    AF=0.19;AC=1;NS=150;AN=2;EAS_AF=0.14;EUR_AF=0.39;AFR_AF=0.09;AMR_AF=0.31;SAS_AF=0.2;VT=SNP;DP=21845     GT:PS   1|0:48005
chr10   106057  .       T       C       .       PASS    AF=0.13;AC=1;NS=150;AN=2;EAS_AF=0.02;EUR_AF=0.0;AFR_AF=0.32;AMR_AF=0.19;SAS_AF=0.01;VT=SNP;DP=22559     GT:PS   0|1:48005

```

Example of truth lines

```
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  sample
chr10   100554  .       C       CT,CTT  30      .       .       GT:AD   1|2:0,1,1
chr10   100588  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   102385  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   102636  .       T       C       30      .       .       GT:AD   1|1:0,2
chr10   102757  .       T       C       30      .       .       GT:AD   0|1:1,1
chr10   104815  .       GCC     G,GC    30      .       .       GT:AD   2|1:0,1,1
chr10   105170  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   105365  .       G       T       30      .       .       GT:AD   0|1:1,1
chr10   106057  .       T       C       30      .       .       GT:AD   1|0:1,1
chr10   106110  .       C       G       30      .       .       GT:AD   1|1:0,2
chr10   106261  .       A       G       30      .       .       GT:AD   1|0:1,1
chr10   108612  .       T       G       30      .       .       GT:AD   1|0:1,1
chr10   108646  .       A       AT      30      .       .       GT:AD   0|1:1,1
chr10   108834  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   110840  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   111743  .       A       AT      30      .       .       GT:AD   1|0:1,1
chr10   112016  .       T       A       30      .       .       GT:AD   1|0:1,1
chr10   112262  .       C       T       30      .       .       GT:AD   1|1:0,2
chr10   113006  .       C       T       30      .       .       GT:AD   0|1:1,1
chr10   113031  .       G       A       30      .       .       GT:AD   0|1:1,1
chr10   113136  .       G       A       30      .       .       GT:AD   1|0:1,1
chr10   113359  .       A       T       30      .       .       GT:AD   1|1:0,2
chr10   113583  .       C       G       30      .       .       GT:AD   0|1:1,1
chr10   113997  .       G       C       30      .       .       GT:AD   1|1:0,2

```

# Answer
Assuming your files are as you show, and the genotype field is always the 1st entry of the sample information field (the 10th), then you can just do this directly with a simple awk script:

```
$ awk -F"\t" -vOFS="\t" 'BEGIN{print "#CHROM","POS","ID","REF","ALT","Truth","Test"}{ if(/^[^#]/){var=$1$2$4$5; gt=substr($10,1,3); if(NR==FNR){a[var]=gt; next;} if(var in a && a[var]!=gt){print $1,$2,$3,$4,$5,gt,a[var]}}}' truth.vcf test.vcf
#CHROM  POS     ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0

```

Here's the script split into separate lines for clarity (you can still copy paste directly into a terminal):

```
awk -F"\t" -vOFS="\t" '
    BEGIN{
        print "#CHROM","POS","ID","REF","ALT","Truth","Test"
    }
    {
        if(/^[^#]/){
            var=$1$2$4$5;
            gt=substr($10,1,3);
            if(NR==FNR){
                a[var]=gt;
                next;
            }
            if(var in a && a[var]!=gt){
                print $1,$2,$3,$4,$5,gt,a[var]
            }
        }
    }' truth.vcf test.vcf

```

The basic idea is:

* `awk -F"\t" -vOFS="\t"`: set the input and output field separators to a tab.
* `BEGIN{print "#CHROM ... }`: print a header.
* `if(/^[^#]/){`: skip the header lines of the input vcf files
* `var=$1$2$4$5`: get a "signature" for this variant, by concatenating the chromosome, position, ref and alt fields.
* `gt=substr($10,1,3);` : get the genotype value.
* `if(NR==FNR){`: if this is the 1st file/
* `a[var]=gt; next`: save the genotype for this variant and skip to the next line.
* `if(var in a && a[var]!=gt){`: this part will only be run for the 2nd file because of the 'next' in the previous block. If this variant was in the truth set and if its genotype there isn't the same as in the test set, print the variant and the two genotypes.

---

If you need the error rate, you can extend the script thusly:

```
awk -F"\t" -vOFS="\t" '
BEGIN{
    print "#CHROM","POS","ID","REF","ALT","Truth","Test"
}
{
    if(/^[^#]/){
        var=$1$2$4$5;
        gt=substr($10,1,3);
        if(NR==FNR){
            a[var]=gt;
            next;
        }
        tot++;
        if(var in a && a[var]!=gt){
            bad++;
            print $1,$2,$3,$4,$5,gt,a[var]
        }
    }
}
END{
    printf "Correct: %s of %s (%s%)\n",tot-bad,tot,(tot-bad)*100/tot;
    printf "Wrong: %s of %s (%s%)\n",bad,tot,bad*100/tot;
}' truth.vcf test.vcf

```

On your example dataset, that prints:

```
#CHROM  POS ID  REF ALT Truth   Test
chr10   100588  .   G   A   1|0 0|1
chr10   102385  .   A   G   0|1 1|0
chr10   105170  .   G   A   0|1 1|0
chr10   105365  .   G   T   1|0 0|1
chr10   106057  .   T   C   0|1 1|0
Correct: 0 of 5 (0%)
Wrong: 5 of 5 (100%)

```
> 3 votes

---
Tags: vcf, haplotypes, phasing
---


# Question
Title: Is there a resource where I can donwload E. coli assemblies grouped by their pathotypes?
I will preface this by saying I'm not a microbiologist, so I apologise for any haziness and I'll be happy to expand.

I would like to download the assemblies of different E. coli genomes and group them by their pathotypes (e.g. EPEC, ETEC, EIEC).

I have used enterobase previously to download salmonella genomes broadly by their subspecies, and it does have a set of E. coli genomes. However, they are mainly annotated by their serotype, and not with any broad nomination of pathogenesis. I think I would be able to use the API to grab these serotypes and map them to a type of pathogenicity but that would probably be quite messy.

Does such a database exist? Most I can find are broadly annotated or a are very detailed, but don't allow access to the raw annotation.

Thanks in advance.

# Answer
If ENA archive contains such information, then you can access it programatically through the EBI api.

Also you can access the NCBI database programatically through their eutils.

In both cases you can access annotations, though not sure if they contain what you are looking for in all the cases.
> 1 votes

---
Tags: genome, public-databases, assembly, microbiology
---


# Question
Title: How can I specify an email address in the SLURM snakemake json configuration file?
I am working on a SLURM-based cluster. I am currently using the following slurm json file:

```
{
    "__default__" :
    {
        "A" : "overall",
        "time" : "96:00:00",
        "nodes": 1,
        "ntasks": 1,
        "cpus" : 24,
        "p" : "LONG",
        "mem": 16G,
        "output": "snakemake%A.out",
        "error": "snakemake%A.err"
    }
}

```

along with the following command:

```
snakemake --jobs 200 --use-conda --printshellcmds  --cluster-config slurm.json --cluster "sbatch -A {cluster.A} -p {cluster.p} -t {cluster.time} --output {cluster.output} --error {cluster.error} --ntasks {cluster.ntasks} --cpus-per-task {cluster.cpus} --mem {cluster.mem}"

```

I would like to be able to specify my email and receive an email when a job finishes. I have tried adding a line in the json file with `"email": my@email.com` and an option to the sbatch command `--email {cluster.email}` but that did not work.

Anyone knows how to do it?

# Answer
> 4 votes
As noted by the comments and the SLURM manual (not so much bioinformatics related?)

You need the following terms in your batch/queuing command:

```
--mail-type=ALL
--mail-user=you@email.com

```

I used "ALL" which will send you *all* reports. If you just want the finish message then you would use "END" and "FAIL" like such:

```
--mail-type=END,FAIL

```

Here are the specifics from the manual if you wish to make a different selection:

**--mail-type=**

> Notify user by email when certain event types occur. Valid type values are NONE, BEGIN, END, FAIL, REQUEUE, ALL (equivalent to BEGIN, END, FAIL, REQUEUE, and STAGE\_OUT), STAGE\_OUT (burst buffer stage out and teardown completed), TIME\_LIMIT, TIME\_LIMIT\_90 (reached 90 percent of time limit), TIME\_LIMIT\_80 (reached 80 percent of time limit), TIME\_LIMIT\_50 (reached 50 percent of time limit) and ARRAY\_TASKS (send emails for each array task). Multiple type values may be specified in a comma separated list. The user to be notified is indicated with --mail-user. Unless the ARRAY\_TASKS option is specified, mail notifications on job BEGIN, END and FAIL apply to a job array as a whole rather than generating individual email messages for each task in the job array.

**--mail-user=**

> User to receive email notification of state changes as defined by --mail-type. The default value is the submitting user.

---
Tags: snakemake, slurm
---


# Question
Title: Finding genes especific to microenvironment
I have RNA-seq .bam files for 3 patients, tumour and its **matched derived model**, namely **organoid**, but I don't have any matched normal sample.

Differentially expressed genes between a tumour and its model mainly come from the presence and lack of **microenvironment** in cancers and organoids respectively and response to stimuli (artificial environment in organoid culture vs cancer.

How I can get the differentially expressed genes specific to microenvironment?

Any intuition please?

# Answer
The genes related to the environment can be identified using a linear model where each organoid is compared to the correspondent original tumor. This is assuming that the environment for all the organoids is the same and that the differences between tumors is just from different patient but are from the same organ/region.
> 1 votes

---
Tags: rna-seq, deseq2, cancer
---


# Question
Title: Download program for ΔG for TM insertion
I am trying to calculate the ΔG of TM insertion for many protein sequences.

Has anyone come across the source code or a programmatic version for the von Heijne ΔG for TM insertion? The web server is at http://dgpred.cbr.su.se/index.php?p=instructions

I can't find a link in the papers or on the site.

# Answer
> 2 votes
The application is stored in github (ElofssonLab/dgpred) with a link to the website.

To find it I searched for: dgpred download github

---
Tags: proteins, sequence-analysis
---


# Question
Title: Restriction enzymes with multiple recognition sequences
In a bioinformatics programming course I am currently teaching (using Python, and approached from a mostly CS rather than a biology background) I assigned a project which involves finding sites where restriction enzymes bind to a DNA sequence. The project is heavily influenced by Chapter 9 of the now classic text "Beginning Perl for Bioinformatics" by James Tisdall. I am using the REBASE database, in particular REBASE Format #9. In this format, each restriction enzyme has a recognition sequence which looks like e.g. `C^YCGRG` where `^` gives information about the exact cut-point and upper-case letters such as `Y` and `R` other than `ACGT` are ambiguity codes. The point of the project is to turn these into regular expressions where cleavage information is discarded and the ambiguity codes are expanded into character classes, and then use these regular expressions to search for matches.

All fairly straightforward, but one very observant student noticed that some enzymes had multiple sequences, for example `M.BssHIII` has a total of 5 sequences: `ACGCGT,CCGCGG,RGCGCY,RCCGGY,GCGCGC`. It is easy to see that these correspond to a total of 11 sequences (after the Rs and Ys are expanded), far fewer than the 144 sequences which would match `VSSSSB` (which is the least ambiguous single sequence capable of matching the above).

My question is, what is going on here? Am I safe in simply replacing `,` by `|` (for "or") in the regular expression, or is there something more subtle going on? Only 8 of the nearly 5000 enzymes have such multiple sequences, and neither the Tisdall book nor even the description of the format (both of which seem to assume single sequences) shed any light on the situation.

# Answer
I would honestly just remove the non-homologous cutting site enzymes from the database. Firstly, these are biology students and require a simple picture to grasp the concept of coding via regex pattern matching. Secondly, in a wet-lab nobody uses these enzymes, I was not even aware the "8" existed until your post. In a teaching lab the best motto is "keep it simple".
> 3 votes

---
Tags: phylogenetics, phylogeny, database
---


# Question
Title: How do I tell featureCounts to ignore fragments that overlap with more than one feature?
I am using featureCounts from the subread package to count the number of fragments per gene in an RNA-seq experiment. I feed it a BAM file with paired-end reads aligned to a reference genome. In the help page of the executable, which I have installed using `conda install subread`, there is a section about overlaps, including for example the `-O` option which will assign the reads/fragments to all overlapping features. What is the default behaviour? Assigning to reads/fragments to the feature with largest overlap? not assigning multi-overlapping reads/fragments?

I would appreciate your help on this.

OS: macOS X

featureCounts version: 1.6.3

# Answer
> 1 votes
If a fragment overlaps multiple genomic features then it's unassignable, which is the appropriate way to handle ambiguous mappings when the resulting counts are then given to tools like DESeq2.

---
Tags: software-usage, subread, featurecounts
---


# Question
Title: Retrieving Metabolic networks from WIT Database
How can I get the dataset for all the 43 metabolic networks for this paper: https://www.nature.com/articles/35036627

It looks like the WIT database does not exist anymore.

# Answer
It seems like there are some files related to the databases at the ftp site: http://ftp.mcs.anl.gov/pub/WIT2/. The related files were last updated on 2002.

I found it via searching at google: `site:*.mcs.anl.gov/ WIT`.

Which reports also a page about the database. There seems to be a project to develop WIT3, because there is a mailing list for developers, maybe you can ask them more details if you have problems with the database.
> 2 votes

---
Tags: networks, data-retrieval, public-dataset
---


# Question
Title: Splitting fasta file into smaller files based on header pattern
I have to split this fasta files into smaller files and write them into individual files my files

```
>lcl|CP000522.1_prot_ABO13860.1_1 [locus_tag=A1S_3471] [protein=hypothetical protein] [protein_id=ABO13860.1] [location=1..957] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13850.1_2 [locus_tag=A1S_3461] [protein=DNA replication protein] [protein_id=ABO13850.1] [location=950..1504] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13851.1_3 [locus_tag=A1S_3462] [protein=hypothetical protein] [protein_id=ABO13851.1] [location=complement(2523..3437)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13852.1_4 [locus_tag=A1S_3463] [protein=YPPCP.09C-like protein] [protein_id=ABO13852.1] [location=3538..4788] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13853.1_5 [locus_tag=A1S_3464] [protein=Cro-like protein] [protein_id=ABO13853.1] [location=5039..5629] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13854.1_6 [locus_tag=A1S_3465] [protein=hypothetical protein] [protein_id=ABO13854.1] [location=complement(6340..6906)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13855.1_7 [locus_tag=A1S_3466] [protein=Resolvase] [protein_id=ABO13855.1] [location=complement(7074..7685)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13856.1_8 [locus_tag=A1S_3467] [protein=hypothetical protein] [protein_id=ABO13856.1] [location=complement(8602..9732)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13857.1_9 [locus_tag=A1S_3468] [protein=putative lipoprotein] [protein_id=ABO13857.1] [location=complement(10072..10374)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13858.1_10 [locus_tag=A1S_3469] [protein=Diaminopimelate decarboxylase] [protein_id=ABO13858.1] [location=complement(10367..10723)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13859.1_11 [locus_tag=A1S_3470] [protein=regulatory protein LysR] [protein_id=ABO13859.1] [location=complement(12076..12444)] [gbkey=CDS]

```

The other pattern is

```
>lcl|CP000523.1_prot_ABO13861.1_1 [locus_tag=A1S_3472] [protein=DNA replication protein] [protein_id=ABO13861.1] [location=1..951] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13862.1_2 [locus_tag=A1S_3473] [protein=hypothetical protein] [protein_id=ABO13862.1] [location=3048..4262] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13863.1_3 [locus_tag=A1S_3474] [protein=hypothetical protein] [protein_id=ABO13863.1] [location=4357..5133] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13864.1_4 [locus_tag=A1S_3475] [protein=hypothetical protein] [protein_id=ABO13864.1] [location=6197..8608] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13865.1_5 [locus_tag=A1S_3476] [protein=secretory lipase] [protein_id=ABO13865.1] [location=8705..9403] [gbkey=CDS]

```

So now my idea is how do i parse and write them into individual files such that CP000522 output written to one file and CP000523 written to another file so forth and so on.

So far what i understand is i have to match the pattern after \>lcl so there are other patterns like "`LN997847`" in the file

Not sure how to proceed tried it in R but failed

it can be done with awk and sed which i tried but i can;t define something that parse all header like takes into account `CP` as well as `LN` .

Any suggestion or help would be highly appreciated . here is my file

# Answer
> 9 votes
Here's a simple awk approach:

```
awk '{if(/^>/){split($1,a,"[|.]")}print >> a[2]".fa"}' Protein_FASTA.txt

```

Or, more concisely, just:

```
awk '/^>/{split($1,a,"[|.]")}{print >> a[2]".fa"}' Protein_FASTA.txt

```

When run on the file linked to in your question, that results in the following files:

```
$ ls
AP014650.fa  CP003848.fa  CP007713.fa  CP012005.fa  CP015122.fa  CP017645.fa  CP018422.fa  CP020594.fa  CP023021.fa  CP024577.fa  CP026712.fa  CP027245.fa  CP030108.fa
CP000522.fa  CP003850.fa  CP007714.fa  CP012007.fa  CP015365.fa  CP017647.fa  CP018678.fa  CP020596.fa  CP023023.fa  CP024578.fa  CP026748.fa  CP027529.fa  CP030109.fa
CP000523.fa  CP003887.fa  CP008707.fa  CP012008.fa  CP015366.fa  CP017649.fa  CP018679.fa  CP021322.fa  CP023024.fa  CP025267.fa  CP026749.fa  CP027531.fa  CU459137.fa
CP000864.fa  CP003888.fa  CP008708.fa  CP012953.fa  CP015484.fa  CP017651.fa  CP019218.fa  CP021327.fa  CP023025.fa  CP026126.fa  CP027121.fa  CP027532.fa  CU459138.fa
CP000865.fa  CP003907.fa  CP008709.fa  CP012954.fa  CP015485.fa  CP017653.fa  CP020573.fa  CP021348.fa  CP023027.fa  CP026127.fa  CP027122.fa  CP027608.fa  CU459139.fa
CP001183.fa  CP003908.fa  CP008850.fa  CP012955.fa  CP015486.fa  CP017655.fa  CP020575.fa  CP021783.fa  CP023028.fa  CP026128.fa  CP027124.fa  CP027609.fa  CU459140.fa
CP001922.fa  CP003968.fa  CP008851.fa  CP012956.fa  CP016296.fa  CP017657.fa  CP020576.fa  CP021784.fa  CP023030.fa  CP026129.fa  CP027179.fa  CP027610.fa  JN377410.fa
CP001923.fa  CP004359.fa  CP010398.fa  CP013925.fa  CP016297.fa  CP018144.fa  CP020577.fa  CP021785.fa  CP023032.fa  CP026339.fa  CP027180.fa  CP029570.fa  LN865144.fa
CP001938.fa  CP006769.fa  CP010399.fa  CP014216.fa  CP016299.fa  CP018255.fa  CP020580.fa  CP021786.fa  CP023033.fa  CP026340.fa  CP027181.fa  CP029571.fa  LN997847.fa
CP002523.fa  CP007578.fa  CP010400.fa  CP014217.fa  CP016301.fa  CP018257.fa  CP020585.fa  CP021787.fa  CP023035.fa  CP026705.fa  CP027182.fa  CP029572.fa  LT594096.fa
CP002524.fa  CP007579.fa  CP010780.fa  CP014292.fa  CP016302.fa  CP018333.fa  CP020589.fa  CP022284.fa  CP024125.fa  CP026706.fa  CP027243.fa  CP029573.fa  Protein_FASTA.txt
CP003501.fa  CP007580.fa  CP010782.fa  CP014293.fa  CP017643.fa  CP018334.fa  CP020593.fa  CP022285.fa  CP024419.fa  CP026708.fa  CP027244.fa  CP030107.fa

```

### Explanation

* `if(/^>/){split($1,a,"[|.]")` : if this line starts with a `>`, split the first field on any occurrence of either `|` or `.` and save the results in the array `a`. Since your header lines all start with `>lcl|`, then the string you are looking for and a `.`, this means that the second value in the `a` array will be your target string.
* `print >> a[2]".fa"` : print (append, `>>`) the current line to a file called "whatever the name of this sequence is" (`a[2]`) and `.fa`. This is run for every line in your input file. ***Note that if you run the same command again, you will need to first delete the files created the first time***. If you don't, because I am using the `>>`, you will just append to the existing files.

---
Tags: fasta, awk, shell
---


title: Filter VCF for "CAF=[0.9,0.1]"

text:
# Question
Title: Filter VCF for "CAF=[0.9,0.1]"
Related to this other question but in this case the CAF is defined as "An ordered, comma delimited list of allele frequencies...".

I tried using vcffilterjdk and the command suggested by Pierre

```
java -jar dist/vcffilterjdk.jar -e 'variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf

```

but predictably returned an error

```
java.lang.NumberFormatException: For input string: "0.1]"

```

I don't know javascript so how should I change the command to make it work?

EDIT:

HEADER

```
##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column.  The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF.  This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

```

LINE

```
chr1    51479   rs116400033     T       A       .       PASS    ASP;CAF=[0.8926,0.1074];COMMON=1;G5;G5A;KGPROD;KGPhase1;OTHERKG;RS=116400033;RSPOS=51479;SAO=0;SSR=0;VC=SNV;VLD;VP=0x050000000005170016000100;WGT=1;dbSNPBuildID=132

```

# Answer
Printing a bracketed list of numerical values is not typical for VCF—I'd guess it's not strictly allowed. A list should be printed only with the commas to separate values. Perhaps the solution is as simple as removing the square brackets.

```
cat input.vcf | perl -ne 's/CAF=\[([^\]]+)\]/CAF=$1/; print' > fixed.vcf

```

Replace `cat` with `gunzip -c` if your VCF file is gzip-compressed.

The Perl regular expression is

* looking for the string `CAF=[...]`
* grabbing the characters in `...`
* replacing the entire match with `CAF=...`

Here's a detailed breakdown of the Perl regular expression.

* `s/.../.../;` — substitution syntax; replace the contents of the left block with the contents of the right block
* `[^\]]+` — one or more non-`]` characters; this is the trickiest bit; see this tutorial for more info on the `[]` syntax for character classes and the special `^` symbol for character class negation
* `(...)` — the parentheses capture the matched pattern into the `$1` variable
* `CAF=$1` — replaces everything matched in the first block with `CAF=` followed by the captured characters
> 0 votes

---
Tags: vcf, file-formats
---


# Question
Title: Definition of "seed" in sequence alignment
I would like to know what is meant by "seed" for various sequence aligners. How is it important?

# Answer
> 12 votes
The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.

# Answer
> 6 votes
Devon's answer gives a good, concise definition. But it's also helpful to consider why seed-and-extend is used and what benefits it provides.

Finding approximate string matches requires operations that are expensive computationally. On the other hand, finding exact string matches can be performed using much less expensive operations—that is, it can be done more quickly with less memory. So quickly finding the location of all seed matches and then applying the approximate matching ONLY AT THOSE LOCATIONS (the seed-and-extend strategy) makes sure that less time is wasted on chasing spurious similarities.

So how long should a seed be? Riffing off your comment, a seed 3bp in length would be entirely ineffective since that same 3bp sequence is going to occur ALL OVER THE GENOME IN RANDOM PLACES: genes, pseudogenes, transposable elements, regulatory sequences, you name it. On the other hand, a seed 51bp in length would also be ineffective since almost every seed would contain at least one sequencing error or variant/polymorphism, thus preventing a perfect match.

In practice, mappers designed for aligning 100-300bp Illumina reads typically use seeds in the high teens or low 20s.

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


# Question
Title: How to interpret amino acid representation
What is the difference between C383S and C383A? For C383S, does it mean Cysteine(C) position 383 in the sequence and its next amino acid connect to it is Serine(S)?

# Answer
This means that at position 383 the Cysteine(C) is mutated to a Serine(S). Depending on the organism and gene usually the mutation notation is linked to a specific transcript since the position could vary depending on transcript splicing.
> 4 votes

---
Tags: proteins, variation, nomenclature
---


# Question
Title: How to interpret UniProt allele patterns?
In creating a parser for the UniProt flat file `uniprot_sprot.dat` (directory) in the manual for the comments `POLYMORPHISM` there are listed patterns for alleles, e.g. for `1A11_HUMAN` is

> The following alleles of A-11 are known:
> A\*11:01 (A-11E), A\*11:02 (A-11K), A\*11:03, A\*11:04, A\*11:05 and A\*11:07.

I tried searching for a description of the patterns using the words `allele` and/or `polymorphism` and found nothing.

* How are the allele patterns to be interpreted?
* Is there a search term to be used when looking for the patterns?
* Also, is there an official repository of the alleles?

The answer does not have to explain the specifics of the example given.

---

Example of explanation of other allele pattern from here

This explanation of the patterns doesn't explain the part of the pattern in UniProt like `(A-11E)`.

# Answer
These are cytotoxic T-cell HLA alleles. HLA genotyping is very common and easy to do, so Genbank is the repository.

A\*11:01 has very high frequency in Aborigine populations here . You can explore the population genetics and the past, present and current population genetics distributions per population of HLA at http://www.allelefrequencies.net/
> 4 votes

---
Tags: sequence-alignment, phylogenetics, phylogeny, file-formats, uniprot
---


title: Matching these matrices in R

text:
# Question
Title: Matching these matrices in R
I have two matrices; I want to convert the row names of first matrix to gene symbol from matched ensemble=gene symbol from second matrix

```
> head(mat1[,1:2])
                TCGA-L5-A4OG-11A-12R-A260-31
ENSG00000000003                         1818
ENSG00000000005                            0
ENSG00000000419                         1436
ENSG00000000457                         1175
ENSG00000000460                          242
ENSG00000000938                          536
                TCGA-IC-A6RE-11A-12R-A336-31
ENSG00000000003                         4596
ENSG00000000005                            3
ENSG00000000419                          751
ENSG00000000457                          840
ENSG00000000460                          205
ENSG00000000938                          253
>

> dim(mat1)
[1] 56925    11
>

> head(mat2)
                ensembl              symbol
ENSG00000274572 ENSG00000274572       ZYXP1
ENSG00000159840 ENSG00000159840         ZYX
ENSG00000162378 ENSG00000162378      ZYG11B
ENSG00000232242 ENSG00000232242    ZYG11AP1
ENSG00000203995 ENSG00000203995      ZYG11A
ENSG00000070476 ENSG00000070476        ZXDC
>

> dim(mat2)
[1] 36848     2
>

```

How I can do that in R?

Thank you

# Answer
Since they share the ensembl\_ID column, you can merge them, then assign the symbol column to the rownames, then delete the symbol column.

something like:

```
merged <- merge(mat1, mat2, by = 0)
rownames(merge) <- merged[,'symbol']
merged[,symbol] <- NULL

```

https://stat.ethz.ch/R-manual/R-devel/library/base/html/merge.html
> 4 votes

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


# Question
Title: Preparing PDB file for modeling with Swiss-Model
I recently starting studying bioinformatics and have a question. I need to build a model of protein using Swiss-Model. First, I need to download a protein record from PDB (for example, 2HZI), then I must open PDB file using simple text editor. I know this usually requires biopython, but i don't use it yet. Now I'm trying to study how to build a model.

Swiss-model needs only one chain to build model, and I remove from file all extra information like AUTHOR,REMARK,HETATM etc. But Swiss-model still gives an error: "Sequence String contains illegal characters. Only letters and gaps allowed."

# Answer
As explained in the website the data introduced must be:

> Format must be FASTA, Clustal, plain string, or a valid UniProtKB AC

So you need to paste just the protein sequence. Without seeing the data you paste I would guess that you don't use a valid plain string.

As Jonathan Moore explains below a FASTA format is:

```
>SEQUENCE_1
MTEITAAMVKELRESTGAGMMDCKNALSETNGDFDKAVQLLREKGLGKAAKKADRLAAEG
LVSVKVSDDFTIAAMRPSYLSYEDLDMTFVENEYKALVAELEKENEERRRLKDPNKPEHK
IPQFASRKQLSDAILKEAEEKIKEELKAQGKPEKIWDNIIPGKMNSFIADNSQLDSKLTL
MGQFYVMDDKKTVEQVIAEKEKEFGGKIKIVEFICFEVGEGLEKKTEDFAAEVAAQL

```

(breaking the sequence in several lines is optional but usually done fore easy reading)
> 3 votes

---
Tags: pdb, proteomics
---


# Question
Title: Pipeline for extracting gene from multiple genomes for use in HyPhy selection analyses?
I have been trying to obtain some preliminary data from HyPhy selection analyses to inform a larger project. I have obtained a number of assembled mammalian genomes from NCBI with the initial goal of extracting from each of them their corresponding sequence for a specific gene. I wanted to use HyPhy to determine if there are signatures of positive selection in any of the mammalian lineages for that specific gene. However, I have been struggling to make headway with this for some time and am realizing that I may be out of my depth.

I was hoping that people could suggest, or at least point me in the direction of, a pipeline/guide/best practices for how one would go about achieving what I am trying to do. I would like to start from scratch and make sure I am performing my analyses in a proper way, using advice from other users.

As I mentioned, I currently have a collection of assembled genomes from NCBI which I would like to pull specific genes from. My understanding is that I will then need to extract the exons (I imagine there is a way to do this without having to extract the entire gene first), align the sequences from each species while preserving reading frames, and then use that alignment as input in HyPhy's codon models. It is very possible that I am missing out some fundamental steps in the process.

If anybody could offer direction/advice, I would be very appreciative.

# Answer
Okay well, because you're using HYPHY, you're going to need a gene tree of every gene you intend to analyze. This requires you to understand the homology relationships between the genomes at hand. Depending on your project, you might already have a list of homologs you're trying to analyze, or you can start from scratch.

If starting from scratch, you will need the proteomes of the species you're interested in (that is, the translation of the coding sequences of each gene, or to put it another way, the concatenated exons). Genome releases should typically have a proteome available for download, along with its corresponding coding sequence.

Next, you need some way to get the homolog groups for the genes in your genomes. You can use a pipeline like OrthoFinder, or iterative tree building and trimming, to infer homolog clusters (OrthoFinder calls them orthogroups but they're the same thing - homologous genes for the species of interest including both single copy orthologs and paralogs).

* For this step, if you already have a set of genes and you know the copy numbers, or know that they're single copy orthologs, then it is more easy - you can just pull the sequence from each proteome and go from there.

Once you have your homologs, you need to align the amino acids using a tool like MAFFT, and then align the coding sequences to the amino acid alignment to attain a coding sequence alignment. For this you can use a tool like pal2nal, or check out `pxaa2cdn` from the `phyx` package written by some of my collaborators.

* If the sequences aren't too divergent, and intron-exon structure is conserved, you could possibly also include alignments of the introns to provide more phylogenetic information. If you do use this approach however, you'd need to align the exons and introns separately. Also, the introns would not be included in the selection analysis.

Next, I would typically reinfer the gene tree of each gene using nucleotide models, which depending on the time, level of selective constraint, etc. can give more signal for the phylogenetic relationships. For this you can use a tree inference program like RAxML. With this tree and your coding sequence alignment, you'll be ready to run analyses in HYPHY (aBS-REL, MEME, FUBAR, etc., depending on your requirements).

* Note that if you want, you can also do this in a Bayesian framework to integrate over uncertainty in the model parameters, tree topology, etc. Several Bayesian phylogenetic inference tools (e.g. MrBayes, PhyloBayes), allow you to infer putatively positively selected sites at the same time as tree inference.
> 2 votes

---
Tags: sequence-alignment, genome, codon, exon, positive-selection
---


title: Attractor Landscape Analysis

text:
# Question
Title: Attractor Landscape Analysis
I have come across a modeling toolbox, ATLANTIS, which is able to determine cell fates in silico based on the input models provided. This MATLAB-based toolbox is built on a method called "Attractor Landscape Analysis." I tried to find a paper dedicated to explaining this method but I did not have any luck so far. Can any of you folks please explain it to me, especially the term "attractor"?

# Answer
From Wikipedia: In the mathematical field of dynamical systems, an attractor is a set of numerical values **toward which a system tends to evolve**, for a wide variety of starting conditions of the system.

In biology you can think of an attractor as a state towards the system is moving to. A good example is the process of differentiation where stem cells compromise to more specialized cells, that could be thought as the attractors of the process.
> 1 votes

---
Tags: networks, algorithms, pathway, literature-search
---


# Question
Title: Affymetrix tags (same ID's) present in different places of the genome
could you explain me how to interpret such situation?

Same tag aligned in several distint places, near different genes... Is the annotation able to distuinguish such 'same' tags, during DGE analysis? Why tags arent unique?

# Answer
> 2 votes
> Is the annotation able to distinguish such 'same' tags, during DGE analysis?

No, you'll end up discarding such probes (assuming you have a reason to actually use microarrays still).

> Why aren't tags unique?

Platforms such as HG-U133 date back to 2002 (meaning they were designed a few years before that), when the human reference genome (and the reference transcriptome) were not nearly as "finished" as they are today. Probes that may have made sense given what we knew almost 20 years ago may no longer makes sense now given what we know. Further, microarrays are rather conservative, in that they tend to primarily gain probes over time and not have them removed (lest companies lose customers who are relying on given probes).

# Answer
> 1 votes
I always assumed that these tags, such as 239963\_at you show, mean that the sequence is not specific of a genome region. This means it cannot be assigned to any gene for DGE.

It could be due to changes on the reference genome or due to a faulty design of the tag sequence. However once in the microarray template they cannot be changed without introducing more problems, that's why doing RNA-seq prevent this.

Sorry I don't have any reference for this (but if someone wants to add this edit the answer).

---
Tags: microarray
---


title: mapping heteryzygous kmers on a genome

text:
# Question
Title: mapping heteryzygous kmers on a genome
I have a set (a couple of millions) of kmer pairs of known length (usually 21) that I know that are heterozygous in the middle nucleotide. For instance:

```
AAAAAAAAAANAAAAAAAAAA

```

Where N has always two states. I would like to map these nucleotides on a fragmented genome assembly that is presumably haploid. I bet that this problem is resolved, but I have hard times to make a practical solution, preferably within python environment. I considered following

1. **Using an existing mapper**

I tried to use python binding of minimap2, thinking that I just need to tell it that I search for exact matches only, but I did not manage to set it up. I set up smaller minimisers but it always returned no hits for kmers of length 21 (I tried to tweak parameters `k`, `w` and `min_cnt`, but to be honest I don't have a good understanding of how minimizers work). Is there a mapper that would be easily set up to do these searches?

2. **Hash table**

The second thing I tried was to build a hash table of all k - 1 / 2 mers of the genome to map the left and right sides of the kmers. Where sequences were mapped to lists of corresponding genomic positions, but this have blown my memory very soon (~50m of kmers were indexed using ~8GB of memory; i.e. this approach does not seems to be viable on a genomic scale).

I am also suspecting that there will be an elegant solution using a suffix array, but I have not cracked that one yet.

# Answer
I don't know if this is the best answer, but one quick-n-dirty approach I commonly use involves tweaking `bwa mem` to report only perfect matches. So in your case, the idea is to make two copies of each 21-mer, replace the N with the appropriate allele, and then map the resulting sequences to your bwa-indexed assembly like so.

```
bwa mem -k 21 -T 21 -a -c 5000 assembly.fasta your-21mers.fasta

```

The `-k 21` says to use a minimum seed length of 21, which forces exact matches. The `-T 21` requires a minimum score of 21, which also enforces an exact match. The `-a` parameter reports all matches, since a "best" match doesn't make sense in this situation. The `-c` parameter limits how many matches are reported, which may need to be adjusted depending on how repetitive the 21-mer is.

**UPDATE** If you don't want to introduce a dependency like bwapy, it's fairly straightforward to implement this in pure Python.

```
with tempfile.TemporaryFile() as kmerfile:
    for kmer in kmers:
        print('>{k}\n{k}'.format(k=kmer), file=kmerfile)
    kmerfile.flush()             # Often works without these two lines
    os.fsync(kmerfile.fileno())  # But best to be safe
    bwacmd = ['bwa', 'mem', '-k', ksize, '-t', ksize, '-a', '-c', '5000', refrfile, kmerfile.name]
    bwaproc = subprocess.Popen(bwacmd, stdout=subprocess.PIPE, universal_newlines=True)
    for line in bwaproc.stdout:
        # process SAM records line by line

```

Also, one can call bwa mem via a C API, but the mechanics of this are left as an exercise for the reader.
> 5 votes

# Answer
Indeed, with suffix arrays, you can map a kmer in ~0.0007s. The way to go is to build a suffix array (for instance using PySAIS package) and then do a convoluted binary search of the four possible kmers. A full implementation of such mapper can be found here.

To build suffix array

```
from Bio import SeqIO
from Bio.Seq import Seq
from PySAIS import sais
from bisect import bisect_right

# load the genome
kmer_genome_file = 'my_genome.fasta'
ffile = SeqIO.parse(kmer_genome_file, "fasta")
sequences = []
scf_names = []
for seq_record in ffile:
    scf_names.append(seq_record.name)
    sequences.append(str(seq_record.seq).upper())

ffile.close()

# paste all scaffolds into a single string and build a suffix array
genome = '$'.join(sequences)
sa = sais(genome)

```

To transform genome coordinate to scaffold/position coordinates we define a function `evaluated2assembly_position`

```
def evaluated2assembly_position(eval_pos):
    scf_index = bisect_right(scf_sizes, eval_pos)
    eval_pos - scf_sizes[scf_index - 1]
    if scf_index > 0:
        return((scf_names[scf_index], eval_pos - scf_sizes[scf_index - 1]))
    else:
        return((scf_names[scf_index], eval_pos))

```

and finally the function for a convoluted binary search (assuming all kmers being 21 nt long and having the unknown nucleotide in the middle, on 11th nt). Note that you need to map to both strands of DNA, which is equivalent of mapping kmer and it's reverse complementary sequence

```
def searchKmer(self, kmer):
    reverse_complementary_kmer = str(Seq(kmer).reverse_complement())
    alignments = []
    for strand in ['+', '-']:
        if strand == '-':
            kmer = reverse_complementary_kmer
        L_kmer = kmer[0:10]
        l = 0
        r = len(sa) - 1
        while l <= r:
            m = (l + r) // 2
            eval_pos = sa[m]
            L_genome_kmer = genome[eval_pos:(eval_pos + 10)]
            if L_kmer == L_genome_kmer:
                r_L_converged = r
                for N in ['A', 'C', 'G', 'T']:
                    kmer = kmer[0:10] + N + kmer[11:21]
                    while l <= r:
                        m = (l + r) // 2
                        eval_pos = sa[m]
                        genome_kmer = genome[eval_pos:(eval_pos + 21)]
                        if kmer == genome_kmer:
                            m_hit = m
                            # following ms
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m += 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            # previous ms
                            m = m_hit - 1
                            eval_pos = sa[m]
                            genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            while kmer == genome_kmer:
                                scf, pos = evaluated2assembly_position(eval_pos)
                                alignments.append([scf, pos, strand, N])
                                m -= 1
                                eval_pos = sa[m]
                                genome_kmer = genome[eval_pos:(eval_pos + 21)]
                            break
                        elif genome_kmer < kmer:
                            l = m + 1
                        else:
                            r = m - 1
                    r = r_L_converged # recover the r after converging on the left sub-kmer
                break
            elif L_genome_kmer < L_kmer:
                l = m + 1
            else:
                r = m - 1
    return(alignments)

```

Taking one of the fist kmers in the scaffold one to test it

```
searchKmer('ACGGCCACCCNGTGTCGTTGT')
Out[23]: [['scaf000001', 28, 'G']]

```

And a softmasked kmer that we can expect more than once

```
kmer = 'ccgaaccacgNtggaggcctg'.upper()

searchKmer(kmer)
Out[31]:
[['CCGAACCACGATGGAGGCCTG', '1_Tps_b3v07_scaf000306', 5196, 'A'],
 ['CCGAACCACGCTGGAGGCCTG', '1_Tps_b3v07_scaf008224', 15781, 'C'],
 ['CCGAACCACGGTGGAGGCCTG', '1_Tps_b3v07_scaf001152', 71261, 'G'],
 ['CCGAACCACGTTGGAGGCCTG', '1_Tps_b3v07_scaf002780', 1721, 'T']]

```

So to map all your kmers (assuming a file with a kmer per line)

```
kmer_file_name = 'my_kmer_file.txt'
with open(kmer_file_name, 'r') as kmer_file:
    kmers = [kmer.rstrip() for kmer in kmer_file]

mapping_list = [searchKmer(kmer) for kmer in kmers]

```

Stuff can be probably done in a more elegant way, but this is a working solution that scales to any genome size.
> 4 votes

---
Tags: genome, algorithms, k-mer
---


# Question
Title: Differential expression analysis when nested effects
I have 3 tumour samples from 3 patients from an experiment, I also downloaded 10 normal samples from TCGA. My design is like this

```
> mycols
               condition batch
normal1           normal     1
normal2           normal     1
normal3           normal     1
normal4           normal     1
normal5           normal     1
normal6           normal     1
normal7           normal     1
normal8           normal     1
normal9           normal     1
normal10          normal     1
OESO_036_a_RNA    cancer     2
OESO_013_a_RNA    cancer     2
OESO_005_a_RNA    cancer     2
>

```

But DESeq2 returns error when trying for differential expression analysis

```
> dds <- DESeqDataSetFromMatrix(countData=counts, colData=mycols, design=~batch+condition)
Error in checkFullRank(modelMatrix) :
  the model matrix is not full rank, so the model cannot be fit as specified.
  One or more variables or interaction terms in the design formula are linear
  combinations of the others and must be removed.

  Please read the vignette section 'Model matrix not full rank':

  vignette('DESeq2')
In addition: Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors

```

Then how I could extract deferentially expressed genes between tumour and normal samples while they are coming from two different batches?

Any help please?

# Answer
llrs is right. "Matrix Model not full rank" means that two of your conditions are confounded. You can't get your normals from one place and your experimentals from another

And no, there generally is no clever way out of this. For all you know, **any** differences you see between your normals and your cancers are 100% due to being sequenced somewhere else, and have nothing to do with real biology.

If you want to run this through DESeq, drop batch effect. But tell whoever you are reporting the results to that batch is hopelessly confounded with cancer state, and you can't tell which differences are due to what.
> 3 votes

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


# Question
Title: Is there a safe catch-all adapter sequence for trimming?
I would like to trim/mark adapters using trimmomatic or picard MarkIlluminaAdapters from a series of Illumina Paired-End read fastqs. The fastq files may have been done using different kits or different generations/versions of illumina chips.

Without knowing ahead of time how the libraries were prepared (or which kit was used), Is there a safe "catch all" adapter sequence (or set of such sequences) that can be passed to a trimming tool and that would cover *all* Illumina kits? In other words, what are we expected to know about the fastq to ensure correct trimming?

*edit:* Picard defines a few built-in adapters which cover many standard kits, but it's not immediately obvious if you need to pass a common one (built-in sequences defined here).

# Answer
You're best off just using `fastp` or `Trim Galore!`, both of which will determine the adapter sequence for you. Trim Galore! uses a built-in list of known sequences for this, whereas fastp uses read overlap.
> 2 votes

# Answer
Alternatively, I really like using bbduk which is part of the BBMap suite.

I've processed every nascent sequencing dataset that has been published, and found a lot of quirky errors with older datasets using TrimGalore. bbduk is a little more fine-tuneable relative to cutadapt/trimmomatic/trimGalore (built on top of cutadapt)/fastp and the run time and memory usage is significantly less.

When you clone/install the repository, there is a directory of adapter fa files that you can choose from for trimming in the resources/ subdirectory which you can specify using the "ref" argument (i.e. ref=adapters.fa). The "adapters.fa" file is an excellent catchall for adapters, but you will likely over-trim a little. Alternatively, if you know what adapters were used in your dataset, you can choose one of the other pre-configured fa files. There are also a number of stats/marking/reporting/trimmed read save options which you can see here:

https://github.com/BioInfoTools/BBMap/blob/master/sh/bbduk.sh
> 2 votes

---
Tags: illumina, picard, trimming
---


title: Tools for quality trimming at 5 prime?

text:
# Question
Title: Tools for quality trimming at 5 prime?
I'm looking for a tool that is capable to do quality trimming at 5' end and it is configurable. For example, I can choose the quality trheshold, the read length, etc. Any recommendations?

I'm currently using sickle for quality trimming. However, one of my samples have very low quality at the beginning as well and sickle doesn't seem to be able to deal well with it.

# Answer
You can have a look at cutadapt. It is capable of quality trimming for both ends as you can read here.
> 4 votes

# Answer
I never spent too much time on choosing my trimming software, therefore I might have missed some jewels. I use trimmomatic when I need versatility and I am entirely sure that it trims reads on both ends. A more convenient option is skewer, they do perform quality filtering, but I could not find an explicit confirmation that it works on 5' too.
> 3 votes

# Answer
Out of all of the major trimming tools available and widely used (trimmomatic, cutadapt/trimGalore (trimGalore is built on top of cutadapt), fastp), I actually instead prefer bbduk which is part of the bbmap suite:

https://github.com/BioInfoTools/BBMap

https://github.com/BioInfoTools/BBMap/blob/master/sh/bbduk.sh

It's much more "fine tuneable" than the other tools and runs in much less time using much less memory (as long as you set the memory limit as it's written mostly in java!). You can trim both ends or use a sliding window approach similar to what you see in Trimmomatic. Even better you can limit the number of bp in that it will trim from either left/right (5'/3'). It's by far the best all-around trimming tool for any type of sequencing dataset ({RNA, ChIP, PRO, GRO, ATAC, DNA}-seq etc...) in my opinion.
> 2 votes

---
Tags: ngs, illumina, quality-control, trimming
---


# Question
Title: How do I prevent the FeatureHeatmap function from the Seurat package, from sorting my data groups in alphabetical order when plotting data?
How can I prevent a function from sorting my data groups (factors) in alphabetical order without affecting the integrity of the data? I am analysing single cell RNA sequencing data using Seurat 2.3.4. I want to be able to visualise the expression of two genes, Gene1 and Gene2 across cell clusters using the `FeatureHeatmap()` function from Seurat. I have assigned to the metadata slot of my Seurat object the factors or conditions being tested in the experiment. However on the resulting FeatureHeatmap plot, the data are sorted and grouped in alphabetical order to my disfavour. I have attempted a solution given on stackoverflow such as the one below:

```
conditions <- c(rep("WTNaClctrl", 1229),
               rep("WTBlu", 1328),
               rep("ILKOBlu", 349))

conditions <- as.data.frame(conditions)

seurat_object@meta.data$conditions <- factor(conditions$conditions, levels=unique(conditions$conditions))

```

When I use the `FeatureHeatmap()` function, the plot **still** sorts the data groups in alphabetical order:

```
FeatureHeatmap(seurat_object, features.plot = c("Gene1", "Gene2"),
               pt.size = 0.25,
               group.by = "conditions",
               max.exp = 3)

```

This plot gets produced:

I would like groupings not to appear in alphabetical order as such:

As previously suggested StackExchange, I have tried to upgrade my Seurat object to Seurat V3:

```
UpdateSeuratObject(Seurat_object):
Object representation is consistent with the most current Seurat version.

devtools::install_github(repo = 'satijalab/seurat', ref = 'release/3.0') : package ‘curl’ successfully unpacked and MD5 sums checked:

Error: (converted from warning) cannot remove prior installation of package ‘curl’.

```

Thank you in advance for any suggestions.

# Answer
> 1 votes
I had this same issue. I reordered the levels and that worked for most functions, but would switch back with the FeatureHeatmap. Turns out there is a line in the function that changes it to a character.

I ended up making a "new" function and just changed that one line from:

```
data.plot$ident <- as.character(x = object@ident)

```

to:

```
data.plot$ident <- object@ident

```

Once that was fixed, the order that I assigned worked out just fine.

---
Tags: r, scrnaseq, seurat, ggplot2
---


title: How to convert featureCounts to FPKM?

text:
# Question
Title: How to convert featureCounts to FPKM?
I have seen many posts regarding counts to RPKM and TPM. I haven't seen any post for counts to FPKM.

I have RNA-Seq data which is paired-end reads. Extracted the counts using featureCounts for all the samples.

There is a function to convert counts to RPKM: using the gene\_length

```
rpkm <- function(counts, lengths) {
  rate <- counts / lengths
  rate / sum(counts) * 1e6
}

```

I know that RPKM is mainly used for single-end reads data. Do you think I can use the above function for converting counts to FPKM as my data is paired end? \[TCGA data: HTSeq counts are converted to HTSeq FPKM\]

Somewhere I have seen that the same function applies for single end data which will be RPKM and also for paired end data which will be FPKM. Am Ir right or wrong?

If not can anyone show some function or code to convert counts to FPKM please.

thanq

# Answer
> 4 votes
> I have seen many posts regarding counts to RPKM and TPM.

There’s your answer then: FPKM = RPKM. It’s simply a more accurate name.

Speaking of RPKM for paired-end data is discouraged because the reference to “read” in this context lends itself to ambiguity. But mathematically the quantity is the same: we are counting *fragments*, not individual reads (of which each fragment has two, for paired-end data).

But as mentioned in the comments, there are good reasons against using FPKMs.

# Answer
> 2 votes
You can use countToFPKM package. This package provides an easy to use function to convert the read count matrix into FPKM matrix; following the equation in

The `fpkm()` function requires three inputs to return FPKM as numeric matrix normalized by library size and feature length:

* `counts` A numeric matrix of raw feature counts.
* `featureLength` A numeric vector with feature lengths that can be obtained using
biomaRt.
* `meanFragmentLength` A numeric vector with mean fragment lengths,
which can be calculate with Picard using CollectInsertSizeMetrics.

See https://github.com/AAlhendi1707/countToFPKM

---
Tags: rna-seq, rna, gene-expression, fpkm
---


# Question
Title: How to extract gene expression tables from this GEO dataset?
I've downloaded this GSE43013 dataset using GEOquery in R. My understanding is that it contains expression data from liver, kidney, and brain for multiple species. I would like to produce gene expression tables for each species and tissue type. The goal would then be to perform PCA and regression analyses to see if there is a relationship between expression of certain genes and e.g. taxonomic position, diet, etc.

However, I have never worked with gene expression data before so I have very little understanding of the GSE file. I have failed to extract any kind of expression data from the file. Solutions that I have found online do not seem to be working for this dataset (for example producing empty tables where I assume there would normally be data on expression levels). Because of my inexperience I cannot determine why they are not working in this situation.

I would appreciate some help in understanding the format of this file, but most importantly, in figuring out how to pull expression data from it, if that is possible.

EDIT: I still haven't managed to pull data from the GSE file, but looking through the samples associated with the data set (such as this one), I can see that there are files formatted like this:

mus.lv.3

ENSMUSG00000000001 26.3376

ENSMUSG00000000003 0

ENSMUSG00000000028 0.501533

ENSMUSG00000000031 0.116078

I assume that the values in the left column refer to specific genes. Perhaps the values in the right column represent expression intensity? If so, would these values be normalized and comparable between species as they are?

# Answer
> 2 votes
Assuming that you have downloaded processed file, according to your link, it is either plain text (compressed) or GTF. If you would have plain text you probably would have posted a snapshot or sample of it, so I assume you have the GTF. According to format description for GTF-files:

* Format consists of one line per feature, each containing 9 columns of data, plus optional track definition lines.
* The 9 columns (or fields) are tab-separated and all fields contain a value (missing/NA-values are usually denoted with ".").
* Fields are:
    1. **seqname** \- name of the chromosome or scaffold; chromosome names can be given with or without the 'chr' prefix. Important note: the seqname must be one used within Ensembl, i.e. a standard chromosome name or an Ensembl identifier such as a scaffold ID, without any additional content such as species or assembly. See the example GFF output below.
    2. **source** \- name of the program that generated this feature, or the data source (database or project name)
    3. **feature** \- feature type name, e.g. Gene, Variation, Similarity
    4. **start** \- Start position of the feature, with sequence numbering starting at 1.
    5. **end** \- End position of the feature, with sequence numbering starting at 1.
    6. **score** \- A floating point value.
    7. **strand** \- defined as + (forward) or - (reverse).
    8. **frame** \- One of '0', '1' or '2'. '0' indicates that the first base of the feature is the first base of a codon, '1' that the second base is the first base of a codon, and so on..
    9. **attribute** \- A semicolon-separated list of tag-value pairs, providing additional information about each feature.

You can find more information on it here and here.

IF you don´t have the file and just have the query output, you probably will find the answer here, where extraction of expression data has been discussed in more detail.

---
Tags: r, bioconductor, gene-expression, gse
---


# Question
Title: RNA-seq analysis of mixed viral/host reads with salmon
I did RNA-seq of mammalian cells infected with pox virus. Now, I have read files which contain both host and virus reads. I want to align the reads both to host and viral genome. I was thinking I could concatenate the host and virus genome into one file and run Salmon against this concatenated genome.

However, salmon recommends transcriptome file for assembly which are not available for viruses. Virus genome are available as GenBank or GFF3 format from NCBI. Is there any way I can concatenate these formats into the format that can be used by salmon? Or is there any way around to use virus genome as reference in salmon?

---

Citing the website:

Salmon is a tool for quantifying the expression of transcripts using RNA-seq data. Salmon uses new algorithms (specifically, coupling the concept of quasi-mapping with a two-phase inference procedure) to provide accurate expression estimates very quickly (i.e. wicked-fast) and while using little memory.

Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods.

# Answer
You provide Salmon with a transcriptome fasta... so merging the human transcriptome with the pox virus genome fasta file should work. You don't mention that in your post but the viral genome should exist as a fasta. You can extract the sequences from genbank or maybe you can find our virus here: https://www.ebi.ac.uk/genomes/virus.html

You can find the coding transcriptome at the link below in the **cdna** folder. The non-coding transcriptome is in the **ncrna**. You should probably merge them:

ftp://ftp.ensembl.org/pub/release-95/fasta/homo\_sapiens/
> 2 votes

---
Tags: rna-seq, covid-19, salmon, virus
---


# Question
Title: Finding new drug-protein targets/interactions using a known one as base
Given a known interaction of a molecule with some proteases (having X-rays images too) is it possible to predict interactions with other proteins either using their 3D structure or not?

If anyone has already tried such techniques I would like to hear a little brief on that subject. Which are the possible techniques, computational power needed and so on.

# Answer
> 1 votes
If you have information on which part of the protein are interacting, and its sequence, you can search for homologous or similar sequences in other proteins. This would give you a list of potential targets and probably a similarity score which you can also use. Then you would have something to run simulations on.

And I think that with 3D structure you should be able to run simulations *in silico* in this case (eventhough it is probably more common to look for the ligands to the protein, not proteins to the ligand).

*In silico* simulations are computationally intensive, but I´ve also seen some being run on a laptop, so you don´t necessarily need access to HPC-clusters, or supercomputer.

---
Tags: protein-structure, interactions, drugs
---


# Question
Title: Where is an up to date miRNA database and what happened to miRBase?
miRBase 21 was published June 26, 2014 and was still in its growth phase. Why is it not being updated anymore or the project declared officially dead? ENSEMBL also uses miRBase as a starting point (http://www.ensembl.org/info/genome/genebuild/ncrna.html).

What do people use now as a good reference for miRNAs especially for less well annotated/sequenced species?

# Answer
> 3 votes
miRbase is a database of miRNA sequences, their genomic locations and the evidence for their existence. It is not a database of their targets, so I am assuming that this is the information you are asking for.

Despite not having been updated since 2014, as far as I am aware miRBase is still the go to place for this data. For example, courses at the European Bioinformatics Institute were still teaching the use of miRBase in April this year (2017).

RNACentral is repository for all non-coding RNAs, but as it is a synthesis of other databases, including miRBase, I suspect data their will still be largely the same as miRBase.

# Answer
> 6 votes
If you are still interested, last year miRBase generated new updates. Currently, according to ftp site the last release is 22.1. So, it is not a dead project and for more specific information you should reference the miRBase blog. In order to get the best set of miRNA annotations, you have to define your target species or even if exists some related annotated species and then, retrieve the current annotation along different miRNA databases, but most importantly you have to be really clear about the definition of miRNA families, because it is clearly different in miRBase and also in RFAM.

# Answer
> 3 votes
I mostly use TargetScan predictions. But here are few curated databases:

* miRGate

> miRGate is a curated database of human, mouse and rat miRNAs/mRNAs targets. It is designed to analyze miRNA and gene isoforms lists under a common and consistent space of annotations. Including all existing 3 UTR and the entirely known miRNAs. All Havana biotypes and ENCODE principal isoforms for the three organisms are also included. Publication
* miR2Disease :

> a manually curated database, aims at providing a comprehensive resource of miRNA deregulation in various human diseases. Publication
* PASmiR

> PASmiR, a literature-curated and web-accessible database, was thus developed to provide a solid platform for the collection, standardization, and searching of these miRNA-stress regulation data in plants. Publication

# Answer
> 3 votes
For miRNA precursor predictions in species in Ensembl but not in miRBase, you can have a go to MapMi (http://www.ebi.ac.uk/enright-srv/MapMi/) .

You can easily run it locally for other species or if you have smallRNA sequencing done on something new.

Let me know if you run into difficulties or have ideas to improve it.

Ensembl itself has its own miRNA predictions.

# Answer
> 2 votes
I would guess DIANA-Tarbase ( evidences, https://www.ncbi.nlm.nih.gov/pubmed/25416803 ) and TargetScan (predictions, https://www.ncbi.nlm.nih.gov/pubmed/26267216 ) are good starting points.

Especially for the less annotated species, TargetScan might be helpful.

---
Tags: sequence-annotation, database
---


# Question
Title: Batch alignment of inconsistently named Fastq files
I have numerous gzip-compressed paired-end Fastq files with ChIP-seq data that I would like to align with bowtie2. I confirmed that bowtie2 takes `.gz` files, although I do see this warning when I run for two samples:

> Warning: gzbuffer added in zlib v1.2.3.5. Unable to change buffer size from default of 8192.

My files are named like so.

```
10_S10_L002_R1_001.fastq.gz  12_S12_L002_R1_001.fastq.gz  2_S2_L001_R1_001.fastq.gz  4_S4_L001_R1_001.fastq-003.gz  6_S6_L001_R1_001.fastq-002.gz  8_S8_L002_R1_001.fastq.gz
10_S10_L002_R2_001.fastq.gz  12_S12_L002_R2_001.fastq.gz  2_S2_L001_R2_001.fastq.gz  4_S4_L001_R2_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz  8_S8_L002_R2_001.fastq.gz
11_S11_L002_R1_001.fastq.gz  1_S1_L001_R1_001.fastq.gz    3_S3_L001_R1_001.fastq.gz  5_S5_L001_R1_001.fastq-002.gz  7_S7_L002_R1_001.fastq.gz      9_S9_L002_R1_001.fastq.gz
11_S11_L002_R2_001.fastq.gz  1_S1_L001_R2_001.fastq.gz    3_S3_L001_R2_001.fastq.gz  5_S5_L001_R2_001.fastq-003.gz  7_S7_L002_R2_001.fastq.gz      9_S9_L002_R2_001.fastq.gz

```

The bash commands I tried are below. When I run this test command:

```
for i in $(ls *.fastq*.gz | grep -v 6 | rev | cut -c 16- | rev | uniq); do
    echo $i;
done

```

...I get this output.

```
12_S12_L002_R1_001.fastq.gz  12_S12_L002_R2_001.fastq.gz

```

When I run this:

```
for i in $(ls *.fastq*.gz | rev | cut -c 16- | rev | uniq); do
    echo $i
    bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 -1 ${i}R1_001.fastq.gz -2 ${i}R2_001.fastq.gz | samtools view -bS - >${i%}.bam
done

```

I do see one BAM for my paired end files. But as I see the second set of files which are named:

```
6_S6_L001_R1_001.fastq-002.gz  6_S6_L001_R2_001.fastq-003.gz

```

if I run the same script over them, I would get an error as the loop can't read the paired end file.

So the tedious way is I can make a separate folder for these sets of data with `fastq-002.gz` and `fastq-003.gz` and can run the script by tweaking the `cut` parameter.

But I would like to know what is the programmatic way of resolving the issue.

# Answer
A very simple approach, but one which will work with the names you show, is to remove any occurrence of `00N` (where `N` is any number) from right before the `.gz`. To do so, edit your script to:

```
for i in $(ls *.fastq*.gz | sed 's/00[0-9]\.gz/.gz/' | rev | cut -c 16- | rev | uniq); do
  bowtie2 -q -p20 -x /run/media/punit/data1/BowtieIndex/hg38 \
    -1 ${i}R1_001.fastq.gz \
    -2 ${i}R2_001.fastq.gz |
        samtools view -bS - >${i%}.bam"
done

```
> 3 votes

---
Tags: read-mapping, fastq, shell
---


# Question
Title: Storing FindAllMarkers results in Seurat object
I am currently working on multiple datasets where each is managed by a separate `Seurat` object. For each I'd like to also compute the marker genes using the `FindAllMarkers` function. Where would it be a good practice to save the markers returned from the `FindAllMarkers` function? I tried adding an extra slot to the `Seurat` object with no success.

I saw the `calc.params` slot and though I could add an item to that list. However, it doesn't seem like the right place as that slot is a list of the parameters used for the different functions executed on the data.

Currently I have some wrapping object that has a slot for the `Seurat` object and a slot for the results from the `FindAllMarkers`. Is there any better place to put it? Somewhere within the `Seurat` object?

# Answer
You can stash anything you like in the `misc` slot (present in both v2 and v3 Seurat objects). I often use it for storing marker information to help with organization.
> 3 votes

---
Tags: seurat
---


# Question
Title: Downloading Fasta for each chromosome of a complete genome
Whatever I am googling I am not able to find fasta file of each chromosome for hs37d5.fa. hs37d5.fa whole fasta file however is here

ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2\_reference\_assembly\_sequence/hs37d5.fa.gz

Do you know where I can get fasta file for each chromosome for this genome?

# Answer
> 2 votes
Personally, I would download the entire file and split it using something like `fastaexplode` from the `exonerate` package. On an Ubuntu machine, you can do that with:

```
wget ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz
gunzip hs37d5.fa.gz
sudo apt install exonerate
fastaexplode hs37d5.fa

```

Alternatively, you could do this on the fly while downloading:

```
wget -O - ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/technical/reference/phase2_reference_assembly_sequence/hs37d5.fa.gz |
  gunzip | awk '/>/{name=$1; sub(/>/,"",name)} {print >> name".fa"}'

```

# Answer
> 2 votes
It is common that reference genomes will be published only as a complete genome, without separate files for each chromosome. Fortunately, a google search for "split fasta by chromosome" will give you many options for downloading a complete genome Fasta file and splitting it by chromosome. The top answer on this biostar thread looks promising.

```
pip install pyfaidx
pyfaidx -x sequences.fa

```

Alternatively, I typically use a homebrew Python script. If there are unlocalized/unplaced scaffolds in the file, this will create a bunch of files with ugly unexpected names. These are easy to ignore or delete.

```
# usage: ./split.py < genome.fasta
from Bio import SeqIO
import sys
reader = SeqIO.parse(sys.stdin, 'fasta')
for record in reader:
    outfile = 'chr{chrid}.fasta'.format(chrid=record.id)
    with open(outfile, 'w') as outstream:
        count = SeqIO.write(record, outstream, 'fasta')

```

With a little bit of searching, you will find a wealth of helpful information on blogs, forums, and library documentation.

---
Tags: fasta, reference-genome
---


# Question
Title: Aligned base strand in pysam `pileupcolumn`
Is it possible that on the same pileupcolumn object (i.e a specific aligned base) I find the same base but aligned to different strands? For example, can I find a coverage of 50 'A's that align to strand '+', and on the same time 50 'A's that align to strand '-' within the same index? I expect that if 'A' is aligned to strand '+', I will find 'T' aligned to strand '-'.

# Answer
> 1 votes
As @terdon mentions, it's certainly possible to have an `A` aligned on both strands. If a read has a mismatch at the location of interest and that mismatch is the complement of the reference base, then you are likely to see `A` in both strands. However, it is very unlikely you will see the same base *at high coverage* in both strands. For example, if your total coverage is 55, and you see 50 `A`s on the + strand and 50 `A`s on the - strand, it's likely that the pileup object is reporting the reverse complement sequence of the reads that align to the - strand for your convenience.

---
Tags: python, pysam
---