# Question
Title: Mapping Reads to Known Gene Paralogs with Long Read Technology
I have some sequencing data from a captured region that is a known paralog **edited**. For now, I have been mapping the data using standard minimap2 flags for PacBio DNA sequencing:

```
minimap2 -ax map-pb GRCh37.fa fastqfile.fastq > samfile.sam

```

I captured a 10kb region and what I am seeing in IGV is puzzling. The region directly before and directly after the paralogues are covered at about 1000X, but the genes themselves are not covered at all:

The two genes in the middle are the paralogues. I left them off to protect data. This data is also only the sorted data, the duplicate marked data has lower coverage (600X).

### Question 1

I was wondering if there is a different approach to mapping I could take that would output **all** alignments for a read, and then specifically only multi-mapped reads?

### Question 2

Is this process normally accomplished by minimap2, or are other aligners more effective for paralogous mapping?

### Question 3

The interesting thing is, while my average read length was 6kb, I am seeing the reads here are about 4kb. I have another region with very high coverage where the reads spanned both paralogues and mapping was successful.

Question: **Why am I only seeing reads that are about 4kb long when the average read length for this dataset is 6kb?**

Any help on this would be great! Thanks and I hope I have been clear on my question,

Dennis

# Answer
### Question 1

This command will get you all of supplementary alignments for the reads. This isn't exactly what you want though. You want all of the reads that have more than one mapping.

```
samtools view -f 2048 -h myfile_sorted.bam > supp_only.bam

```

This bash script returns a bam file that contains all alignments for reads that have multiple mappings. It works by only selecting reads that have supplementary alignments (2048 flag), or are not the primary alignment (256 flag), then using picard tools to filter the `bam` file. Using picard tools is the fastest way I have seen to filter a `bam` file given a list of reads.

```
#/bin/bash
BAMFILE=/path/to/my_file.bam
OUTFILE=/path/to/multiply_mapped.bam
PICARD_JAR=/path/to/picard.jar
TEMP=multiply_mapped.txt

for filt_int in 256 2048; \
do samtools view -f "$filtint" ${BAMFILE} | \
  grep -v '#' | cut -f1; done | sort | uniq > ${TEMP}; \
java -jar ${PICARD_JAR} FilterSamReads \
  I=${BAMFILE} \
  O=${OUTFILE} \
  READ_LIST_FILE=${TEMP} \
  FILTER=includeReadList
rm ${TEMP}

```

This command took 21 seconds to execute for my bam file with ~300,000 long reads mapped and several repetitive regions.

### Question 2

Filtering the `bam` file should be pretty standard, as addressed above in my response to **Q1**. However, different long-read aligners are more sensitive to distantly-related paralogs. You might consider trying the blasr aligner since it is tuned for PacBio long reads, or try bwa-mem to compare your results with minimap2. Realistically though, most of these aligners will probably behave similarly to one another.

**The fact that there are paralogs does not matter to the aligners.** The aligner just sees these regions as two different places to which a single read can map just about equally well. When there are multiple similar regions to which a read can map, as in the case of paralogs, the aligner will report both alignments. **I do not know of any long-read paralog-aware aligners.**

As to why none of the reads mapped to the paralogs: Perhaps they have low mappability, as in they are repetitive? Also, maybe the capture didn't work properly? It looks like there is serious sequencing bias to the dataset given the uneven coverage.

### Question 3

Are the reads that are mapping actually 4kb, or are the alignments only 4kb long? You may have many 6kb-long sequences that have been clipped at either the 5' or 3' end. If this is not a heavily-curated genome then the assembly may be incorrect around the paralogs, preventing proper mapping.

One suggestion is to take whole-genome long reads from the same organism and to map them. Then check if the mapping to this region is equally strange compared to your capture data. This will give you a better idea of what caused the uneven coverage and seeming lack of coverage around the paralogs.
> 3 votes

# Answer
I have found an alternative solution that has been excellent for this paralog mapping issue, so I thought I'd post what I had worked on:

**Question 1+2**

I was wondering if there is a different approach to mapping I could take that would output all alignments for a read, and then specifically only multi-mapped reads?

Using blasr, pacbio specified some options in their github docs for long, nearly identical reads:

```
blasr input.bam /path/to/reference.fa --hitpolicy all --bam --out alignments.bam --minMatch 8 --maxMatch 15 --nproc 4

```

The link to the docs as to why this is useful is here, although I did tweak the parameters to smaller seed sizes to create more anchors to the reference, which produced more secondary alignments. This is effective because this will output up to 10 secondary alignments per read, which I saw quite frequently for my region of interest.

After this, sorting was accomplished using Picard SortSam, with the following flags:

```
java -jar picard.jar SortSam INPUT=${SAMPLE} OUTPUT=./sorted/${SAMPLE}\
SORT_ORDER=coordinate VALIDATION_STRINGENCY=LENIENT\
CREATE_INDEX=true MAX_RECORDS_IN_RAM=100000

```

The blasr files have much more information than normal sams produced by aligners like minimap2 and bwa mem, so the MAX\_RECORDS\_IN\_RAM=100000 as opposed to 500000 helps to ensure less frequent job failures (this was done on a cluster).

Next, I parsed out the region of interest:

```
samtools view -h input.bam "chrX:1000-2000" > extracted.bam # this is fake

```

Parse out mapq=0 reads:

```
awk '$5==0' extracted.bam > mapq0.txt

```

Get read names to grep for later:

```
awk '{print $1}' mapq0.txt > read_names_mapq0.txt

```

Search for all reads with readname from the original bam to get mapping positions:

```
samtools view infile.bam | grep -f read_names_mapq0.txt > all_reads_mapq0_mapq_high.txt

```

Keep only columns: 1) read\_name, 2) Sam Flag, 3) Chr, 4) Pos, 5) MapQ:

```
awk '{print $1,$2,$3,$4,$5}' OFS='\t' all_reads_mapq0_mapq_high.txt > mapq_read_info.txt

```

Then I wrote a simple python script to sort by read name, count the number of times a read shows up, and see all the positions it maps to. An example of the output:

```
        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       272     19      3395674 0
read1/10158955/ccs       0       16      217682  254
read1/10158955/ccs       272     12      64990527        0
read1/10158955/ccs       256     16      222913  0
read1/10158955/ccs       272     19      21328060        0
read1/10158955/ccs       256     19      3224052 0
read1/10158955/ccs       272     17      34878621        0
read1/10158955/ccs       272     16      84204312        0
read1/10158955/ccs       256     1       244788680       0
read1/10158955/ccs       272     2       222951342       0
read2/10223868/ccs       16      16      217057  254
read2/10223868/ccs       256     20      3078724 0
read2/10223868/ccs       272     6       133521134       0
read2/10223868/ccs       256     7       152515722       0
read2/10223868/ccs       256     21      45258117        0
read2/10223868/ccs       256     1       249160685       0
read2/10223868/ccs       272     16      223631  0
read3/10551872/ccs       272     GL000220.1      144982  0
read3/10551872/ccs       256     18      72831954        0
...

```

read\_counts output:

```
        Read_Name       Flag    Chr     Pos     MAPQ
read1/10158955/ccs       10      10      10      10
read2/10223868/ccs       7       7       7       7
read3/10551872/ccs       9       9       9       9

```

The numbers in each column correspond to the number of supplementary reads seen. A brief explanation of sam flags which I feel is relevant for unpaired reads:

```
Flag=0: read was mapped to + strand
Flag=16: read was mapped to - strand
Flag=256: read was mapped to + as secondary alignment
Flag=272: read was mapped to - as secondary alignment

```
> 0 votes

---
Tags: long-reads, pacbio, mappability, minimap2, paralogs
---


title: scRNA-seq, 10x cellranger pipelines

text:
# Question
Title: scRNA-seq, 10x cellranger pipelines
I'm starting to do sc-rnaseq using 10x cellranger pipelines, and i add TdTomato sequence to mouse reference genome and add an entry in the gtf. My code is But when I makref, it remindered that: Writing genes GTF file into reference folder... Invalid number of columns in GTF line Can anyone tell me the reason?

# Answer
As was suggested by @benn, you will need to add the TdTomato sequence to the end of your fasta formatted genome file. If your TdTomato sequence is in fasta format already, `cat` will do the trick. If not, just go the end of the file, start a new line and add a header line starting with "\>" then immediately the name of your gene, no whitespcaes in between, something like "\>TdTomato". Then add your sequence in a new line.

The format of the GTF file is explained in the image you posted, the trick is to use the exact same name (in column 1) that is used in the genome file you will create at the previous step (no "\>" though). Again, this line, in which the different fields (columns) will be separated by tabs, will be added to the end of your GTF file.
> 2 votes

---
Tags: scrnaseq, fasta, gtf
---


# Question
Title: What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?
What are the pros and cons of the different basecallers in Oxford Nanopore Technology Sequencing?

I am about to start a MinION run on my laptop. What should I consider when choosing my basecaller? Can I let MinION do its sequencing and generate fast5 files and call the bases later? I think I might use Albacore. Metrichor and MinKNOW can also do basecalling. Are there reasons to choose one over another?

# Answer
> 6 votes
First of all - yes, you can generate FAST5 files and basecall later. Basecalling during the sequencing run is useful if you want results more quickly. You can also recall your FAST5 files with multiple basecallers, if desired.

There are several ways to basecall currently:

1. MinKNOW
2. Albacore
3. Guppy

MinKNOW uses an embedded version of Albacore to perform its basecalling. However, the version in MinKNOW occasionally lags behind the released Albacore version.

Guppy is a GPU enhanced basecaller, and has not yet been released by ONT. It is used by default on the GridION. Results are similar, but not identical to Albacore, but Guppy is much faster if you have a supported GPU.

There used to be a cloud version of basecalling in Metrichor, but this is no longer available, and all basecalling must now be performed locally.

I'd suggest looking over Ryan Wick's comparison of the different available base callers, as he has done a comprehensive look (mainly from a qualitative perspective) at the different basecallers.

# Answer
> 7 votes
There are also third party free and open source basecallers that haven't been developed by Oxford Nanopore.

Of particular note is Chiron, which gave the best uncorrected assembly identity among the base callers that Ryan Wick tested. It's slower than Albacore, but appears to be more customisable, and could theoretically be used to model any DNA feature that is detectable in the electrical signal, rather than just the A/C/G/T base sequence.

There's a paper about Chiron on BiorXiv.

# Answer
> 5 votes
It's worth noting two things as of Dec 2018:

1. Albacore is being deprecated (but is still available from the Nanopore developer portal).
2. Guppy is under active development, so Ryan Wick's comparisons may not reflect the current state of things. (The claim is that the current, just-released version of Guppy uses a new "flip-flop" algorithm that improves accuracy over Albacore.)
3. MinKnow contains the current version of the production basecaller. As of the release today, this has switched from Albacore to Guppy. If you want, you can run MinKnow for your basecalling (which will run Guppy), but you will likely need a beefy computer with access to a good amount of storage to do that.

# Answer
> 1 votes
Depending on your local system I'd suggest to run the MinION without local basecalling, and perform basecalling later using Albacore on server/cluster. That will keep storage requirements on your nanopore-computer low. Have albacore output fastq files which will be okay for most typical applications. Keep the fast5 files to repeat basecalling with newer versions of basecallers and to call nucleotide modifications.

---
Tags: nanopore, ngs, minion, genome-sequencing, base-calling
---


# Question
Title: Are there free cloud computing platforms for biology projects?
I want to implement the analysis found in the paper RNA-Seq of Tumor-Educated Platelets Enables Blood-Based Pan-Cancer, Multiclass, and Molecular Pathway Cancer Diagnostics The project id is 281708. It has 285 samples comprising about 300 Gb of read data. The reads are available on sra.

I want to apply a machine-learning pipeline to find important features in these data but I do not have access to computing resources to finish the project.

## Question

Are there free cloud computing platforms that I can use to perform this research?

# Answer
> 1 votes
If you're getting started with a project and need some free computational resources you should look at https://galaxyproject.org/. There are some free servers that you may be able to use. However, Galaxy is designed to have pre-fabricated pipelines and a GUI to avoid having to deal with command-line tools. Therefore any sort of custom analysis will have to go through the galaxy API.

Otherwise, most of the free computational resources that were available several years ago like iPlantCollaborative are no longer available. Running servers costs money, after all, and the funding for such projects often runs out after a few years.

You may consider asking researchers at your local universities to collaborate and provide you with access to their servers.

---
Tags: fastq, genomics, phylogenetics, sratoolkit
---


# Question
Title: Extracting genes from corrplot and adding labels based on high and low corelation
I am doing a gene-gene pairwise correlation, I get a plot and I do see positive and negative correlation

My code:

```
library(corrplot)
library(Hmisc)
gen = read.csv('NEW_RBP/CORPLOT/TF_NORMAL_NEW_40.txt',row.names = 1)
head(gen)
dim(gen)
gendat <- t(gen)
#gendat <- gen
head(gendat)
dim(gendat)

macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
p_mat <- cor_5$P
row_names <- rownames(M)

pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")

corrplot(M, method = "color",is.corr = FALSE,
         tl.srt = 45,diag = TRUE,tl.pos='n'
         ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)

dev.off()

```

How do I extract those gene pairs? One way what I have read and done is:

```
library(dplyr)
flat_cor_mat <- function(cor_r, cor_p){
    #This function provides a simple formatting of a correlation matrix
    #into a table with 4 columns containing :
    # Column 1 : row names (variable 1 for the correlation test)
    # Column 2 : column names (variable 2 for the correlation test)
    # Column 3 : the correlation coefficients
    # Column 4 : the p-values of the correlations
    library(tidyr)
    library(tibble)
    cor_r <- rownames_to_column(as.data.frame(cor_r), var = "row")
    cor_r <- gather(cor_r, column, cor, -1)
    cor_p <- rownames_to_column(as.data.frame(cor_p), var = "row")
    cor_p <- gather(cor_p, column, p, -1)
    cor_p_matrix <- left_join(cor_r, cor_p, by = c("row", "column"))
    cor_p_matrix
}

#cor_3 <- rcorr(as.matrix(gendat))

my_cor_matrix <- flat_cor_mat(cor_5$r, cor_5$P)
head(my_cor_matrix)

dim(my_cor_matrix)

Delete.na <- function(DF, n=0) {
      DF[rowSums(is.na(DF)) <= n,]
    }

TF_NORMAL_CORR = Delete.na(my_cor_matrix)
dim(RBP_DISEASE_40_CORR)

TF_NORMAL_HIGH = TF_NORMAL_CORR[ TF_NORMAL_CORR$cor > 0.9 & RBP_DISEASE_40_CORR$p < 0.00000001, ]

TF_NORMAL_LOW = TF_NORMAL_CORR[TF_NORMAL_CORR$cor < -0.9 & TF_NORMAL_CORR$p < 0.00000001, ]

```

Now even if I do that sort of high cut-off filtering I still get considerably high number of correlation pairs like few thousands if positive correlation is considered.

So few of my issues which I would like to get input or help

* How do I extract the cluster of regions with positive correlation or negative correlation ?
* If I found certain gene pair/pairs how do I label them ? should I make a vector or map or any other way like if say more than 0.8 correlation I would like to put a label on those regions only
* Considering label the gene order is important to be considered, what I have read and found is this how to find the order of the label

> ```
> tempCor <- cor(gendat)
> tempCor <- data.frame(tempCor)
> dim(tempCor)
> names(tempCor) <- c(1:741)
>
> out <- corrplot(t(tempCor),method = "color",is.corr = FALSE,
>                  tl.srt = 45,diag = TRUE,tl.pos='n'
>                  ,order = "hclust",hclust.method = c("complete"),tl.cex = 3,number.cex=1.5)
>
> dimnames(out)
>
> ```

Now with all these how can I extract those patches of high and low correlation pairs and add selected label or highlight them and add label

I m not sure if all these can be done with corrplot! like adding selected label if yes i would be really glad to use that..

My data file data used

"In which groups are you interested in? In the upper left region and the bottom right square? In which format do you want it ? A list for each group? BTW, how many samples do you have? If it is so high, or they are mostly similar or there are few samples."

Its symmetric so would it matter if i use upper left or bottom right? I guess anyone would do good.

Format - a list for each group if i understand if there is one gene with many positive correlation then that would be considered as group?is that you are saying ?

I have 16 samples in this case but here I using gene-gene. I'm bit curious with a naive question: Is it possible to do gene to sample correlation or if i know what genes are positively or negatively correlated can I trace back without have to plot another heatmap of found genes which are positively or negatively correlated to see their expression in sample ?

I have given the data file it might help you better with my context..

# Answer
Here is a solution using the pheatmap library to cluster and visualise the correlation matrix, then extract the groups from the cluster dendrograms:

```
gen <- read.csv(file = "TF_NORMAL_NEW_40.txt", header = TRUE, row.names = 1, sep = ",")
gendat <- t(gen)
library(corrplot)
library(Hmisc)
macolor = colorRampPalette(c("navyblue", "white", "red"))(100)

cor_5 <- rcorr(as.matrix(gendat))
M <- cor_5$r
str(M)
p_mat <- cor_5$P
row_names <- rownames(M)

library(pheatmap)
pdf("TF_NORMAL_NEW_40_gene.pdf",20,15,paper = "a4r")
plot_data_pheatmap <- pheatmap(M, color = rev(macolor), clustering_method = "complete", fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

```

```
tree_cut <- cutree(plot_data_pheatmap$tree_row, h = 20)
tc <- data.frame(tip = names(tree_cut), clust_membership = as.character(unname(tree_cut)))
row.names(tc) <- tc$tip
tc <- tc['clust_membership']
head(tc)
str(tc)
class(tc$clust_membership)
pdf("TF_NORMAL_NEW_40_gene_clustcut.pdf",20,15,paper = "a4r")
pheatmap(M, color = rev(macolor), clustering_method = "complete", annotation_row = tc, annotation_col = tc, fontsize_row = 0.8, fontsize_col = 0.8)
dev.off()

```
> 4 votes

---
Tags: r, visualization, correlation
---


# Question
Title: PCA plot shows big difference but not many differentially expressed genes are found
I got a `PCA` plot of `bulk RNA-seq` experiment that looks the following way:

It was generated by the following code:

```
pcaData <- plotPCA(rld_sva, intgroup=c("Group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed()

```

First `sva` correction was run to correct for batch effects and then `rlog` transformed values were plugged in to `plotPCA` function.

The first issue that catches the eye is that `100%` of the variance is explained by just 2 dimensions. I am not sure what can one say about the data in this case. The second issue is that I get only around `20` differentially expressed genes by using `DESeq2` analysis (`log2FoldChange > 1`, `p_adj < 0.05`). I know that we can not directly state that if there is a large difference on `PCA` there will be present a plenty of differentially expressed genes, but why is it not the case? Simple logic tells me that `pca` shows the difference between the samples in their gene expression, so I would expect seeing a plenty of differentially expressed genes.

# Answer
You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.

On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.
> 7 votes

---
Tags: scrnaseq, deseq2, pca
---


# Question
Title: How to link GDC ids to CCLE cell line names?
Hello bioinformaticians,
I've recently downloaded few RNAseq data files from Genomic Data Commons Data (GDC) portal. These files belong to Broad institutes CCLE project.

Now the problem is that GDC portal do not provide CCLE cell line names:

So my question is how could I link the GDC ids (prefixes of File Names) to CCLE cell line names?

I would imagine because it is an important step in the analysis, there should be some way to do that.

# Answer
The filenames of those bams contain the cell line names. They are constructed:

SAMPLEID.CELL\_LINE\_NAME.REPLICATE.bam

Thus the first 5 cell lines in your screen capture are SK-N-F1, TE-9, CJM, HCI-H1915 and CAL-33.

All the information that is need is also present in the experimental metadata `xml` file attached to each sample record. This is accessible by downloading it from the record on the website, or by querying the GDC metadata API (https://gdc.cancer.gov/developers/gdc-application-programming-interface-api).
> 2 votes

# Answer
I would recommend using UCSC Xena instead of GDC for CCLE data. They don't have everything. For example, BAM files are not available. However, the data that is available is organized in standard text files with consistent identifiers.
> 1 votes

---
Tags: rna-seq, data-download, cell-line
---


# Question
Title: No counts for added gene in cellranger (scRNA-seq)
I have a set of scRNA-seq samples enriched with FACS for cells expressing a specific gene reporter (TdTomato). In particular the gene I want to report has positive counts in the resulting matrix for 97% of cells.

I followed CellRanger documentation to create a custom reference. I added TdTomato on an additional chromosome and I added an entry in the gtf:

```
TdTomato    .   exon    1   1431    .   -   .   gene_id "TdTomato"; transcript_id "TdTomato"; gene_name "TdTomato"

```

The command I'm using is:

```
cellranger-2.1.1/cellranger mkref --genome=refdata-cellranger-mm10-2.1.0.tdtomato --fasta=genome.tdtomato.fa --genes=genes.tdtomato.gtf

```

where `genome.tdtomato.fa` is the original reference genome (mm10, as provided by cellranger) with an additional chromosome at the end with header:

```
>TdTomato

```

and `genes.tdtomato.gtf` the original reference genome gtf with the additional line indicated previously at the end of the file.

I actually tried both on the '+' and on the '-' strand, however I guess this should not affect the read mapping, since TdTomato is on an independent chromosome.

I ran `cellranger count` using the custom reference:

```
s="SAMPLE"
cellranger-2.1.1/cellranger count --id="$s"_tt --transcriptome=refdata-cellranger-mm10-2.1.0.tdtomato --fastqs=original_fastqs_directory/"$s"/fastqs --sample="$s"

```

I am using the default parameters of CellRanger.

However, I found that less than 1% of cells seems to express TdTomato. If I use the '+' strand I get slightly more counts (2.4% of cells have positive TdTomato).

What I checked:

* With IGV I can see that there are reads mapping to TdTomato
* Using samtools I am able to see the same number of reads mapped on TdTom either with the + or - strand
* However, following this documentation I noticed that by executing the CellRanger filtering steps, more reads are filtered in the case of the - strand

This reduces the number of reads used for counts to 44 (2.4% of the original 1807 reads, - strand), 240 (13.3%, + strand) for TdTomato, 71620 (21.2%) for the gene reported by TdTomato. Interestingly, there are fewer reads for TdTom than counts, so there may be some kind of normalization afterwards?

The gene reported by TdTomato is longer, 11394 bp than TdTom, 1431 bp. The total length of the exons for the two genes is 5349 vs TdTom: 1431. Does this influence exponentially the amount of reads mapped to the two genes?

# Answer
The first column of your GTF must correspond exactly to the name of the sequence in your FASTA file. Typically this is the chromosome so a GTF line:

```
chr1    hg38_knownCanonical exon    11106535    11262507    0.000000    .   .   gene_id "gene1"; transcript_id "tx1";

```

This must have a corresponding sequence in the FASTA file:

```
>chr1
ATCG ATCG ATCG ATCG ATCG ATCG...

```

Therefore you need to add an additional sequence to your FASTA file with `>TdTomato`. This must have exactly the same sequence as the construct that you are expecting reads for. The run `cellranger mkref` and `cellranger count` on these FASTA and GTF files. Cellranger also can take a GTF with “CDS” as a feature in addition to “exons”. The STAR alignment algorithm used by cellranger is splice aware so it can handle split reads but you will need to ensure the GTF corresponds to exons of the mRNA/cDNA rather than the full sequence.

10x Genomics is a platform with a high dropout rate, especially for lowly expressed genes. You should bear this in mind for what you expect to observe. It may be worth comparing the mapping rate of your construct to housekeeper genes to see if this is the reason for the low counts. Since whole mRNA molecules are captured on the beads before preparing the library, there should not be a difference in UMI counts between long and short genes.
> 1 votes

---
Tags: scrnaseq, cellranger, tdtomato
---


# Question
Title: Multiple sequence pairwise global alignment against reference sequence
### Background

I have a fasta file with protein sequences. The first sequence in the fasta file is the reference sequence and every subsequent sequence is something I want to align to the reference. The end goal is to identify every unique sequence and be able to compare the regions of variation.

### What I have tried

I'm looking into different tools such as Clustal Omega and MAFFT, but there is no option to select the reference sequence against which all the other protein sequences will be aligned.

### Questions

**Q1:** I want perform one pairwise global alignment of each sequence in the file against the first (reference) sequence. **What software can I use to do this?**

**Q2:** **What specific commands can I use to achieve Q1?** I'm looking for specific instructions on how to complete this, not just a suggestion of a tool, since I'm not able to find anything from a relatively comprehensive search of available tools.

# Answer
### Q1

From this phrase, "The end goal is to identify every unique sequence", I am not 100% sure of what you mean by "every unique sequence". It sounds like maybe you want to group alignments of homologs and visualize them together?

If the proteins are similar enough, you might also benefit from trying BLAST to look at similar regions.

Otherwise MAFFT, Clustal, or MUSCLE are all appropriate tools to generate the alignments that you need. However I also do not know of any software aside from Geneious that has the option to align many proteins to a reference. If you were working with DNA I'm sure the stock suggestion would be to use a read aligner like `bwa mem` or `minimap2` and to use a bam file.

Since I don't know of any software to do exactly what you would like, I wrote a small python script to run MAFFT many times for each sequence in a fasta after the first sequence, the reference.

### Q2 - MAFFT/Python

In any case, here is a bit of python code that performs a MAFFT alignment of all of the sequences in a fasta file against the first sequence in the fasta file, which you have designated as the reference.

Here is the fasta file from Wikipedia. I called it `seqs.fasta`:

```
>reference
MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVSLRDLFDRAVMVSHYIHDLSS
EMFNEFDKRYAQGKGFITMALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHL
VTEVRGMKGAPDAILSRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDED
ARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC*

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGYVLPWGQMSFWGATVITNLFSAIPYIGTNLV
EWIWGGFSVDKATLNRFFAFHFILPFTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLG
LLILILLLLLLALLSPDMLGDPDNHMPADPLNTPLHIKPEWYFLFAYAILRSVPNKLGGVLALFLSIVIL
GLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQMASILYFSIILAFLPIAGX
IENY

```

Here is a script, `mafft_ref.py`. Run it with `python mafft_ref.py seqs.fasta`:

```
from Bio import SeqIO
from Bio.Align.Applications import MafftCommandline
import tempfile

fasta = sys.argv[1]

reference = 0
count = 0
with tempfile.NamedTemporaryFile() as temp:
    for record in SeqIO.parse(fasta, "fasta"):
        if count == 0:
            reference = record
        else:
            SeqIO.write([reference, record], temp.name, "fasta")
            mafft_cline = MafftCommandline(input=temp.name)
            stdout,stderr=mafft_cline()
            print(stdout)
        count += 1

```

Here is the output of `python mafft_ref.py seqs.fasta` printed to the terminal. Use `python mafft.py prots.pep > output.txt` to save the output to `output.txt`.

```
>reference
------------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQ
VSLRDLFDRA--------------VMVSHYIHDLSS-----EMFNEFDKRYAQGKGFITM
ALNSCHTSSLPTPEDKEQAQQTHHEVLMSLILGLLRSWNDPLYHLVTEVRGMKGAPDAIL
SRAIEIEEENKRLLEGMEMIFGQVIPGAKETEPYPVWSGLPSLQTKDEDARYSAFYNLLH
CLRRDSSKIDTYLKLLNCRIIYNNNC
>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
ADQLTEEQIAEFKEAFSLFDKDGDGTITTKEL------------GTVMRSL---GQNPTE
AELQDMINEVDADGNGTIDFPEFLTMMARKMKDTDSEEEIREAFRVFDK---DGNGYISA
A---------------------------------------ELRHVMTN-----------L
GEKLTDEEVD-------EMIREADIDG-------------------DGQVNYEEFVQMM-
-----TAK------------------

>reference
------------MDSKGSSQKGSRLLLLLVVSNLLLCQGVVSTPVCPNGPGNCQVS----
--LRDLF--------------------DRAVM----VSHYIHDLSSEMFNEFDKRYAQGK
GFITMALNSCHTSSL-------PTPEDKEQAQQTHHE-------------VLMSLILGLL
RS--WNDPLYHLVTEVRGMKGAPDAILSRAIEIEEE--------------NKRLLEGMEM
IFGQVI-------------------P-------------------GAKETE-PYPVWSGL
PSLQTKDEDARYSAFYNLLHCLRRDSSKIDTYLKLLNCRIIYNNNC
>gi|5524211|gb|AAD44166.1| cytochrome b [Elephas maximus maximus]
LCLYTHIGRNIYYGSYLYSETWNTGIMLLLITMATAFMGY----VLPWG----QMSFWGA
TVITNLFSAIPYIGTNLVEWIWGGFSVDKATLNRFFAFHFILP-----------------
-FTMVALAGVHLTFLHETGSNNPLGLTSDSDKIPFHPYYTIKDFLGLLILILLLLLLALL
SPDMLGDPDNHM----------PADPLNTPLHIKPEWYFLFAYAILRSVPNK--LGGVLA
LFLSIVILGLMPFLHTSKHRSMMLRPLSQALFWTLTMDLLTLTWIGSQPVEYPYTIIGQM
ASI----------LYFSIILAFLPIAGXIENY--------------

```

References: Some credit of `mafft_ref.py` to a biostars answer found here.
> 2 votes

---
Tags: sequence-alignment
---


# Question
Title: Visualize similar interactions in two homologous networks
## Background

I have two networks with homologous genes from two different species.

To identify homologs I performed a functional comparison. One of the species is not well studied while the other one is the house fly (D. melanogaster). I used ClueGO to predict the functions of the genes in the D. melanogaster network.

I found a function of interest influenced by a set of genes and predicted that the same genes in the other network (thus far uncharacterized) should also perform a similar function.

## Questions

**Q1:** How do I check if there are overlapping patterns between both the networks on Cytoscape?

**Q2:** There are more than 100 homologs in this list and I would like to minimize the list by only considering those homologous genes that are arranged the same as well.

# Answer
> 4 votes
An easy way to visualize this is to make a tab-delimited table that contains edge information for both species and both networks, color the edges depending on the species, and make the thickness and transparency dependent on the interaction strength.

For example, say you have a list of orthologs for the two species and just assign them a number or a name. Each edge has a "source" and a "target". The edge weight is the interaction strength.

(Note do not copy-paste the below table as the mixed tabs/spaces may make the graph topology incorrect). In my toy dataset the data look like this:

```
source  target  variable    value
gene1   gene2   speciesA    3
gene1   gene3   speciesA    42
gene2   gene3   speciesA    3
gene1   gene2   speciesB    9
gene2   gene3   speciesB    5
gene2   gene4   speciesB    41
gene3   gene4   speciesB    32

```

Import your network to cytoscape as a network from a file:

Then go to the style tab:

... and make sure edge is selected.

Then you can adjust the transparency, edge color, and edge stroke settings, like so using edge stroke color:

:

... and transparency panels as examples.

**Don't forget to save your style when it is complete!**

When you have adjusted these settings you will have something that looks like this:

# Answer
> 0 votes
I think what you're looking for, as I suggested in this old question (quite similar to this one, I believe), is the GASOLINE app for Cytoscape that allow you to perform an analysis that will tell you which parts of the two networks are conserved and which genes belong to the overlapping subnetworks.

---
Tags: visualization, networks
---


# Question
Title: Filter Trinity transcriptome based on RNASeq reads
I have recently generated a genome-guided transcriptome with Trinity, and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is with the goal of trying to reduce the initial dataset down to something a bit more manageable (I have about 300k transcripts covering 250Mb in total, but would prefer about 1/10 of that number).

I see that Trinity has a workflow using Bowtie2 for evaluating the quality of assembled transcripts, but Bowtie2 is getting a bit old in the mapping world now. Are any of the super-fast transcript mappers (e.g. Salmon or Kallisto) appropriate to use for transcript filtering?

# Answer
> 6 votes
I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:

1. Run BUSCO in short/transcript mode on the Trinity-generated sequences
2. Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
3. Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
4. Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
5. Sub-select the transcripts based on this threshold, using the Salmon results table

## Step 1 -- BUSCO

```
python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10

```

## Step 2 -- Salmon

```
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A

```

## Step 3 -- Results merge (using R)

```
#!/usr/bin/Rscript
busco.table <-
    read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
               header=FALSE, stringsAsFactors = FALSE,
               col.names=c("Busco id","Status","Sequence","Score","Length"));
busco.table <- subset(busco.table, Status != "Missing");

count.table <- read.delim("quant/reads-all_quant_ISR/quant.sf",
                          stringsAsFactors=FALSE);

busco.expr.df <- merge(busco.table, count.table,
                       by.x="Sequence", by.y="Name", all.x=TRUE);

```

## Step 4 -- Count thresholding (visual exploration)

```
pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
options(scipen=10);
busco.numReads <- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                              function(x){
                                  exp(mean(log(x[x>0])))}), decreasing=TRUE);
par(mfrow=c(1,2));
plot(x=seq(0,1,length.out=length(busco.numReads)),
     busco.numReads, log="y", ylab="Number of Mapped Reads",
     yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
     xleft=0.9, xright=par("usr")[2],
     col="#80808020", border=NA);
plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
     busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
invisible(dev.off());

```

Looks like a value of 50 will be suitable.

## Step 5 -- Sub-select the transcripts

```
cat(file="HighCount_Transcripts_Num50.txt",
    subset(count.table, NumReads >= 50)$Name, sep="\n");

```

Then filter using a suitable FASTA selection utility:

```
pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta

```

# Answer
> 3 votes
I would recommend Transrate. It generates a score for each transcript based on many metrics (read their paper here) and, using that score, gives a predicted cut-off to produce the optimal transcript set.

# Answer
> 0 votes
DETONATE offers a pipeline for filtering in this way as well. But the aligner used by default is Bowtie2...

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


# Question
Title: Multithread fastq processing with kseq.h in C++11?
## Background

I am using the wonderful `kseq.h` as my fastq/fasta parser for a C++ project. My actual implementation uses something like this, where I am counting kmers of each sequence and performing a bunch of other processing on every read pair.

Thanks to @user1782818's code my program is very fast, but I want it to be faster. The only way I know how to do this is with multithreading. However, I have never worked with multithreading in C++.

For each read pair I am accessing a vector or map of structs. Within each struct I am incrementing counters or pushing values onto vectors. Essentially I am doing something like this:

```
//compiled with g++ 7.3 using command g++ this_prog.cpp -o this_prog -lz
// execute with ./this_prog my_file.fastq.gz
#include <iostream>//cout
#include <map>
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
KSEQ_INIT(gzFile, gzread)

int main(int argc, char **argv)
{
  //initialize prefix
  std::map<std::string, int> prefix_counter;
  std::string prefix;
  prefix.resize(2);

  gzFile fp;
  kseq_t *seq;
  int n = 0, slen = 0, qlen = 0;
  fp = gzopen(argv[1], "r");
  seq = kseq_init(fp);
  while (kseq_read(seq) >= 0){
    ++n;
    //COUNT PREFIXES
    std::copy(seq->seq.s, seq->seq.s + 2, begin(prefix));
    if (prefix_counter.count(prefix) == 0){
      prefix_counter[prefix] = 1;
    } else prefix_counter[prefix]++;
  }
  std::cout << "Total n. of sequences: " << n << '\n';
  kseq_destroy(seq);
  gzclose(fp);

  //print counts
  for( const auto& sm_pair : prefix_counter ){
    std::cout << sm_pair.first << ": " << sm_pair.second << '\n';
  }
  return 0;
}

```

Example output from the above program is something like this for a 100,000-read test fastq file.

```
Total n. of sequences: 100000
AA: 8561
AC: 5178
AG: 7634
AT: 8523
CA: 7526

```

## Question

**Q1:** Is it possible to perform a task like the one shown above using multithreading and `kseq.h` to read in a fastq file and process the reads in chunks? I only want to use standard libraries. C++11 is OK - I know that multithreading support was added in C++11.

**Q2:** Can you please provide a C++ example of how to convert the above or similar program to use multithreading using `kseq.h`?

# Answer
You might take a look at my FQFeeder library (https://github.com/rob-p/FQFeeder), which pairs the kseq parser with a multi-producer / multi-consumer queue to allow efficient parallel processing of the parsed reads. The library contains an example program to demonstrate how to use it with multiple threads. It's worth noting that kseq is already *very* fast, so multi-threading only really makes sense if worker threads are doing non-trivial work on each read so that you can then make the parsing and processing concurrent (have the processing keep pace with kseq).

I'll note that my library relies on a header-only MPMC implementation and so technically isn't "only using the standard library", but you're already relying on `kseq.h`, so omitting third-party libraries must not be a hard constraint.
> 7 votes

---
Tags: fastq, c++, multithreading
---


# Question
Title: How can I extract the longest N isoforms per gene from a fasta file?
The question of how one can extract the longest isoform per gene from a multi-fasta file has been answered expertly in this thread: How can longest isoforms (per gene) be extracted from a FASTA file?

I wonder: how can the python script "filter.py" by Devon Ryan, that is marked as best answer, be modified so that it extracts the N longest isoform sequences per gene instead of only the single longest isoform?

Any help on this to a python novice is very much appreciated.

Here is an example of how the data looks (two genes with 5 isoformes each, and I need the 3 longest isoforms for each gene):

```
>Animal1|gene000001_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000001_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000002_1
TTLQLLNELVLLLHYTTLLLHTKKTKESTINNYFKQKEQK
>Animal1|gene000002_2
CYYYTTQHYYYIQKKQKNLQSTITSNKKNKNKHAK
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

```

As output, I would like to have a fasta file like the following:

```
>Animal1|gene000001_5
EILCYSIPPGVGLLFTPVNFGMLRSKLIKVRILDIINFNCSLLICLPPYLFQIDVGEDKD
DPVVEIVLYLACSI
>Animal1|gene000001_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000001_4
HHAVNFQPMSEATGIYLNSTEIRNSLLQYSSWSWFTFHTCKLWNVTF
>Animal1|gene000002_3
ARDMLDYTTTTKRVSATTTLHNTTTTYKKNKRIYNQQLLQTKRTKISMQNKNLNQINYKK
LQETTEKHAK
>Animal1|gene000002_4
LCSAVARQKTKSRKYDATVDTCFSNIKININVLSFSLIVNPKYSKVSYQYGTNDAFRKFF
ITSRSLR
>Animal1|gene000002_5
ETLEYFGLTINEKESTLILILMLLKQVSTVASYFLDLVFCLATALHS

```

# Answer
> 0 votes
Here a python solution:

```
#!/usr/bin/env python
from __future__ import print_function
from Bio import SeqIO
from argparse import ArgumentParser

def main():
    args = get_args()
    last_gene = None
    transcripts = []

    for record in SeqIO.parse(args.input, "fasta"):
        gene = record.id.split("_")[0]

        if last_gene is None:
            last_gene = gene

        if last_gene == gene:
            transcripts.append(record)
        else:
            print_longest_transcripts(transcripts, args.number)
            last_gene = gene
            transcripts = [record]

    print_longest_transcripts(transcripts, args.number)

def print_longest_transcripts(transcripts, number):
    longest = sorted(transcripts, key=lambda x: len(x.seq), reverse=True)

    for record in longest[:number]:
        print(">"+record.id, record.seq, sep="\n")

def get_args():
    parser = ArgumentParser(description="")
    parser.add_argument("-i", "--input", help="input fasta", required=True)
    parser.add_argument("-n", "--number", help="number of longest transcript per gene", type=int, required=True)
    return parser.parse_args()

if __name__ == '__main__':
    main()

```

Run it like this:

```
python longest.py -i input.fasta -n 3

```

# Answer
> 1 votes
I know you are asking for a python solution. But maybe you are interested in other ways.

You can use `samtools faidx` to index your fasta file. The resulting file contains the information about the length of every sequence. We can sort by length and extract the first n transcripts per gene.

**1. Index fasta**

```
samtools faidx input.fa

```

**2. Extract IDs of the n longest transcripts per gene**

```
sed 's/_.\+//g' input.fa.fai | sort -u |
    parallel -k 'grep {} -m 3 <(sort -k2,2 -r --sort g input.fa.fai) |
        cut -f1' > longestID.txt

```

The `sed` command replaces everything from `_` onwards in each line, leaving behind only the gene name. `sort -u` is used to get a list of distinct gene names. For each of theses gene names we `grep` for the first 3 matches in the descending sorted by length index file. As we just want the sequence ids, we select only the first column with `cut`. Adopt the `-m` parameter if you want more or fewer transcripts.

**3. Get sequences**

```
samtools faidx input.fa -r longestID.txt > output.fa

```

---
Tags: rna-seq, filtering, isoform
---


# Question
Title: Why do NEBNext indexing primers have sequence between the p5 oligo and index?
In a previous post I asked **Why do NEB adapters have non-complementary sequence?**

Since then, I realized that there is some other sequence in the p5 indexing primer, as well as in the p7 indexing primer.

Here is a diagram of the NEBNext protocol. The parts that I am confused about are in step 5.

1. My first question is about the p5 index. Why is there a sequence `5'GATCTACAC 3'` between the p5 flowcell oligo and the index?

2. My second question is about both indices. Why is there no gap between the p7 index and the flowcell annealing oligo?

For reference, I fetched these indexing primer sequences from NEB on page 18 and 19 of NEBNext Multiplex Oligos for Illumina instruction manual.

# Answer
Question 1: The additional sequences is needed because that is complimentary to the P5 sequence anchored to the flow cell. It is also the site for priming the Index 2 read on a MiSeq. The additional sequence is also needed for the read 1 primer in the cartridge to anneal to the correct place for the molecules.

Question 2: The "GT" sequence you're referring to is necessary for the Index 1 and Read 2 primer to anneal to the correct place.
> 5 votes

---
Tags: phylogenetics, illumina, adapter, pcr, library
---


# Question
Title: python does not quit when the input file is too big
I have a script that extracts the set of unique kmer pairs from a kmer dump file (& then using them downstream to guess ploidy).

The script seems to work correctly (it's reporting a meaningful set of kmer pairs), but since we have added a parallelization it does not quit when the job is done for really big data sets. The problem is that I never managed to reproduce the behavior with a small reproducible example and it takes around 400G of memory and a day of computation to get to the problem. My imagination is quite limited, I have no idea how should I debug such a problem.

I don't really get it; the script actually prints the last line

> \*_families\_2.tsv and \*_coverages\_2.tsv files saved.

, but then it runs further (& forever). I thought it could be related to longer time to free the memory as discussed here, but it does not seem like the case because it does not quit even after a really long time (I tried to just leave it for a day or so and when i kill it the memory is freed in like 20s).

Does anybody an idea where the problem could be? Or how to figure out?

The reported problem on GitHub.

# Answer
You've hit the RAM bottleneck. Its very common in de novo assemblies and resolved using template assemblies. I don't know this program but suspect it is performing a pairwise analysis for every piece of data in a genomics data set.

The solution is to either:

> 1. find a bigger machine, e.g. a big multi-core machine and make sure the program distributes across every core.
> 2. parallelise the calculation (an excellent way to overcome RAM problems), or,
> 3. find an approximation to the calculation that is less computationally demanding, templates are a good way to circumvent this.
> 2 votes

---
Tags: python, k-mer
---


# Question
Title: How to export web NCBI tBLASTn results in table format with many queries?
### Context

I'm an MSc student working on writing up my thesis (back home now) from my laptop and, therefore, unfortunately don't have access to a workstation/server capable of doing the tBLASTn search that I wanted to do. As a result I have been trying to use as many online resources as possible.

### Question

Is there a way to export the tBLASTn results from the NCBI website to a simple text file? I want the output to be something similar to:

```
FASTA_Name Max_score Total_score Query_cover E_value Identity Accession Gene_ID

```

### Question Details & requirements

Essentially, I have a yeast database with gene expression information (where the gene names are named in the *S. cerevisiae* convention eg YMR105C or YDL204W) and I'm trying to relate each of those expression profiles to their respective homologous gene in my species of yeast.

I've exported the systemic file names (and their respective protein residue sequences) for *S. cerevisiae*, turned that into a big fasta file (~6,000 sequences) and then run a tBLASTn against my yeast genome to give a single Report page where you select each submitted FASTA sequence from the "Results for:" dropdown and it brings up the alignments for that specified sequence as normal. To manually curate this would take forever and I was wondering if there's a way to generate a table containing all the results where the E-value is below a certain threshold? Ideally it'd be something along the lines of:

"FASTA Name" "Max score" "Total score" "Query cover" "E value" "Identity" "Accession" "Gene\*"

\*where "Gene" is what you get when you get the title of the page if you were to click that specific "Accession" hyperlink... appreciate that might be beyond the scope of this question though.

Is this possible (or indeed is there a completely different strategy) or should do I need to download the relevant databases and run something along the lines of this?

Thanks for all your help in advance!

# Answer
> 5 votes
It sounds like your easiest option is to just install blast on your local machine and set up the databases there. It really will only take a few minutes to blast 6,000 sentences against a yeast genome on a modern laptop. Any solution interfacing with web results or parsing the output will be painfully convoluted. There are executables for blast that will work on Mac, Linux, or Windows.

Some beginning instructions for installing it on your laptop can be found here: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE\_TYPE=BlastDocs&DOC\_TYPE=Download . Make sure to pass the --num\_threads option to speed up the search.

---
Tags: public-databases, phylogenetics, blastp
---


# Question
Title: Why does this human bam file only have one copy of each chromosome?
As we know that in human DNA sequence, one copy of chromosome comes from mother's DNA and another copy comes from father's DNA so as to form two copies of each chromosome in human DNA. So, if we extract exome sequence from its DNA then each exome must also have two copies of each chromosome (one from mother and another from father). But in either fastq file or bam file in WES data, I always find only one copy of each chromosome.

**Q1:** Where is the other gene copy in the sequence or have I have missed something?

**Q2:** How can I check the ploidy and find if there are two copies of each chromosome in WES data? How can I do this if my WES data are mapped to the reference and stored in a bam file?

# Answer
> 16 votes
The maternal and paternal copies of a chromosome are called **haplotypes**. Many metazoans (animals) are diploid and have maternal and paternal chromosome contribution during sexual reproduction, not just humans as your question states.

## Response to Q1

Your question, in other words, is: Why do `.bam` files not differentiate between haplotypes?

Your question gets to a more fundamental kernel of how most people do "genomics" today. When most people assemble a genome or create a reference genome, they actually just produce a fasta file that contains a single sequence for each chromosome. This, of course, is not biologically accurate as there are two very distinct sequences per chromosome. Most genomes only report a chimeric consensus between the two and call this a reference. These reference genomes are haploid, or haplotype-collapsed.

**This is where your question comes in:** You have mapped transcripts that biologically are derived from two different haplotypes to a reference genome that only contains one consensus sequence that (inadequately) represents both haplotypes. **As a result, at a single locus in the bam file there will be mapped reads from both haplotypes.** If your reference fasta file contained the exact sequence of both haplotypes, if you mapped reads to the reference, and if you only look at primary alignments, the reads will mostly map to their haplotype-of-origin.

This gets into another topic called **phasing**, wherein the order and orientation of most of the polymorphisms unique to each haplotype are determined using sequencing data. There are some problems with this as it relies on properly detecting variant sites. Software like GATK and others can find single nucleotide polymorphisms (SNPs) if there is good sequencing coverage, however detecting insertions and deletions is much more difficult. This gives a very SNP-biased view of the haplotype differences in any given genome. After finding variant sites, phasing software like hapcut2 determines which variants fall on which haplotype and outputs blocks of variants that are thought to belong to the same haplotype.

Phasing alone is not enough to accurately reconstruct the exact sequence of both haplotypes due to the inability to detect all variants with read mapping. The gold standard of the future is **diploid de novo genome assembly**, in which each haplotype is assembled independently. This is an active area of research for people who develop genome assemblers and is tightly linked to the advances in PacBio and Oxford Nanopore long reads. This paper on trio canu is a good start for learning about one successful dioploid assembly technique.

## Response to Q2

If you want to check ploidy of all of the chromosomes then you need at least 10X coverage of whole-genome shotgun data and a tool like smudgeplot and genomescope.

If you are trying to check local duplications or whole-chromosome duplications for something like a cancer sample you would also need whole-genome shotgun data. WES data does not give reliable information on ploidy since the read count at a given locus is dependent on the transcription level, not the actual quantity of homologous chromosomal DNA for the region in question.

If you are trying to (use a bam file to) look for evidence of chromosomal duplication or duplication of a locus in cancer samples you would look for an increase in whole genome shotgun coverage at that locus compared to a known normal locus. For example if all chromosomes but one have a coverage of 22, and one chromosome has an average coverage of 33, this is evidence of a trisomy. This logic can be applied to smaller regions as well (barring repetitive regions, paralogs, et cetera).

---
Tags: sam, ngs, fastq, exome
---


# Question
Title: Truncating branch length values of Phylogenetic tree with biopython
I have been using biopython 1.72 to display my phylogenetic tree files.

1. Using the function `'Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)'` to display branch lengths as well on tree, the tree displays the branch lengths as including all numbers after decimal. While I want the numbers to be truncated to two decimal places. Following is a python method that truncates upto `"the number required"` decimal places.

def truncate(n, decimals=0): multiplier = 10 \** decimals return int(n * multiplier) / multiplier truncate (12.567, 2) #Using function 12.56 #Output Can I use such method with biopython for branch lengths?
2. `'Phylo.draw'` function opens the window which is very narrow and if the tree file is larger, all the branches are overlapped and it gets difficult to read what is written. Does anyone know any better display method?

Following is the example code, I have been using:

```
from Bio import Phylo
from Bio.Phylo.TreeConstruction import *
from Bio import AlignIO

aln = AlignIO.read(open('example.phy'), 'phylip')

calculator = DistanceCalculator()
dm = calculator.get_distance(aln)

constructor = DistanceTreeConstructor()
njtree = constructor.nj(dm)
starting_tree = njtree
scorer = ParsimonyScorer()
searcher = NNITreeSearcher(scorer)
constructor = ParsimonyTreeConstructor(searcher, starting_tree)
pars_tree = constructor.build_tree(aln)
Phylo.draw(pars_tree, branch_labels=lambda c: c.branch_length)

```

I have an msa file as input and it gives a tree, which looks like:

```
                                                         _______ Human
                                                 _______|
                                                |       |_______ Chimpanzee
                                         _______|
                                        |       |        _______ Dog
                                 _______|       |_______|
                                |       |               |_______ Cow
                         _______|       |
                        |       |       |_______________________ Elephant
                 _______|       |
                |       |       |_______________________________ Mouse
                |       |
         _______|       |_______________________________________ Platypus
        |       |
        |       |        _______________________________________ Anole_lizard
        |       |_______|
  ______|               |        _______________________________ Chicken
 |      |               |_______|
 |      |                       |_______________________________ Zebra_finch
 |      |
 |      |_______________________________________________________ Xenopus
_|
 |       _______________________________________________________ Zebrafish
 |______|
 |      |_______________________________________________________ Fugu
 |
 |______________________________________________________________ Amphioxus

```

# Answer
> 1 votes
Apologies, try the following,

```
allclades = list(pars_tree.find_clades(order='level'))
for i, blength in enumerate(allclades):
        printf(blength=%.2f, blength)

```

"allclades" is a dictionary rather than a simple array.

If this fails I would examine the allclades data structure either using `pprint` (directly) or using `printStruct` via the following code,

```
>>> d = [{'a1':1, 'a2':2, 'a3':3}, [1,2,3], [{'b1':1, 'b2':2}, {'c1':1}], 'd1', 'd2', 'd3']
>>> printStruct(d)

def printStruct(struc, indent=0):
  if isinstance(struc, dict):
    print '  '*indent+'{'
    for key,val in struc.iteritems():
      if isinstance(val, (dict, list, tuple)):
        print '  '*(indent+1) + str(key) + '=> '
        printStruct(val, indent+2)
      else:
        print '  '*(indent+1) + str(key) + '=> ' + str(val)
    print '  '*indent+'}'
  elif isinstance(struc, list):
    print '  '*indent + '['
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ']'
  elif isinstance(struc, tuple):
    print '  '*indent + '('
    for item in struc:
      printStruct(item, indent+1)
    print '  '*indent + ')'
  else: print '  '*indent + str(struc)

```

Once you know the structure you unravel it accordingly to get to the branch lengths and then issue a printf command to round them to 2 dp.

Goodluck

# Answer
> 0 votes
The universal coding way to do this is the `printf` command regardless of language (sprintf is used language depending), for Python this is,

```
blength = 0.22342
printf("Taxa A, blength=%.2f", blength)

```

---

**Output**

Taxa A, 0.22

---

The following may be useful to your context:

```
import numpy
allclades = list(pars_tree.find_clades(order='level'))
for blength in enumerate(allclades):
        printf(blength=%.2f, blength)

```

Goodluck

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


# Question
Title: Question on nanopore sequencing data process pipeline (cDNA-PCR)
I recently started doing the analysis on nanopore sequencing data. As I was searching for some help on pre-processing of the data, I found your nice setup pipeline created here: https://www.protocols.io/view/stranded-mapping-from-long-reads-n8cdhsw.

As I was trying to follow the instruction, I ran into a question on part where it uses lastal to divide the input file as strand specific output.

Part that uses the following script below:

```
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \   awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \   ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \   awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

```

In here, I was wondering what the ONT\_barcodes\_adapters.fa file actually looks like. I wanted to find the information in other places and could not find it elsewhere.

If you could help me with this, I would really appreciate it!

# Answer
> 1 votes
I recommend that you look at the most recent version of the protocol, which includes example adapter and barcode files:

https://www.protocols.io/view/stranded-mapping-from-long-reads-syheft6?step=2

\[check in the protocol to see if there's a more recent version\]

These files contains the public sequences that ONT have released for their sequencing kits (as briefly mentioned in the "Before Starting" section of the protocol). The sequences are available on the ONT community, but I have converted them into FASTA files for use in my protocol.

The adapter and barcode files were combined into a single FASTA file in the previous version of my protocol. For the PCR cDNA barcoding, it looked something like this:

```
>ONT_VNP
ACTTGCCTGTCGCTCTATCTTCTTTTT
>ONT_SSP
TTTCTGTTGGTGCTGATATTGCTGCCATTACGGCCGGG
>BC01
AAGAAAGTTGTCGGTGTCTTTGTG
>BC02
TCGATTCCGTTTGTAGTCGTCTGT
>BC03
GAGTCTTGTGTCCCAGTTACCAGG
>BC04
TTCGGATTCTATCGTGTTTCCCTA
>BC05
CTTGTCCAGGGTTTGTGTAACCTT
>BC06
TTCTCGCAAAGGCAGAAAGTAGTC
>BC07
GTGTTACCGTGGGAATGAATCCTT
>BC08
TTCAGGGAACAAACCAAGTTACGT
>BC09
AACTAGGCACAGCGAGTCTTGGTT
>BC10
AAGCGTTGAAACCTTTGTCCTCTC
>BC11
GTTTCATCTATCGGAGGGAATGGA
>BC12
CAGGTAGAAAGAAGCAGAATCGGA
>RB12A
GTTGAGTTACAAAGCACCGATCAG

```

I separate out adapters and barcodes in the newer versions of the protocol because I was getting too many incorrect hits.

---
Tags: nanopore
---


title: Does Prokka do six-frame translations?

text:
# Question
Title: Does Prokka do six-frame translations?
Does Prokka do six-frame genome annotations? The corresponding publication does not mention it.

> Seemann, Torsten. 2014. “Prokka: Rapid Prokaryotic Genome Annotation.” Bioinformatics 30 (14): 2068–69.

# Answer
Yes, Prokka will consider all six-frames, the gene prediction is actually done with Prodigal, paper: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2848648/

github: https://github.com/hyattpd/Prodigal
> 7 votes

---
Tags: sequence-annotation
---


# Question
Title: How to normalise scRNASeq data for differential expression analysis
I wish to perform differential expression analysis for cluster-specific gene expression in single-cell data (with a tool such as MAST or SCDE).

I have data for 3 biological replicates. I performed analysis on 10x data with cellranger and have filtered out cells with low UMIs. I have corrected for batch effects and identified clusters by performing graph-based clustering on the UMAP projections.

As I understand it, the components from batch effect correction algorithms are *only* intended to be used for the purpose of plotting and clustering. They are *not* to be used for analysis of individual gene expression or differential expression. This has been stated unequivocally by the original authors of MNN and CCA batch correction algorithms.

I have derived a UMI count matrix for each sample which corrects for amplification bias. Differences between the cells and samples still persist in the raw data. These are obvious from the mean UMI counts to the UMAP of uncorrected data. The number of UMIs differs considerably between samples (batches) and individual cells within a sample (the total UMI per cell). The amount of data from each of the cells varies considerably. The read-depth between samples is also different. Therefore, I need a simple approach to correct for differences in read-depth, total UMI or sensitivity between cells.

Is it possible to normalise UMI counts between cells? Should UMI count data be treated differently to reads per transcript in bulk RNA-Seq? Would it be more appropriate to do so at the level of individual cells or for batches across samples?

# Answer
I would suggest using a likelihood ratio test for differential expression using logistic regression with batch as a latent variable. In Seurat you can do:

`markers <- FindAllMarkers(object, test.use = 'LR', latent.vars = 'batch')`

(change "object" and "batch" accordingly)

See https://www.biorxiv.org/content/early/2018/02/14/258566 from Lior Pachter's group for some more info on the LR test.
> 3 votes

---
Tags: scrnaseq, differential-expression, single-cell, normalization, 10x-genomics
---


# Question
Title: Separate boxplots for multiple violin plot
I am using the following function from seurat package to generate multiple violon plots and I am interested in adding box plots to them but it doesn't work when I have plotted different data at once. Is there a way to solve it ?

For example, this works:

```
library(Seurat)
VlnPlot(object = pbmc_small, features.plot = 'PC1') +
  geom_boxplot()

```

But this will simply lead into an empty box on top of my plots:

```
VlnPlot(object = pbmc_small, features.plot = c('PC1', 'PC2')) +
  geom_boxplot()

```

# Answer
> 2 votes
This happens because the violin plots are combined using `cowplot::plot_grid` before being returned by VlnPlot. You can prevent the plots from being combined by setting `combine=FALSE`, then modify each one by adding a boxplot, then combine the modified plots using `Seurat::CombinePlots`.

However, the combine argument is currently broken in VlnPlot. I have pushed a fix but it's not on the public branch yet. VlnPlot is just a wrapper around ExIPlot (expression by identity plot) in Seurat v3 so right now this will work:

```
library(Seurat)
library(ggplot2)

plots <- ExIPlot(object = pbmc_small,
                 features = c('PC1', 'PC2'),
                 pt.size = 1,
                 combine = FALSE)

for(i in 1:length(plots)) {
  plots[[i]] <- plots[[i]] + geom_boxplot() + theme(legend.position = 'none')
}
CombinePlots(plots)

```

Note this answer is for Seurat v3

# Answer
> 2 votes
The `vioplot` package comes built in with boxplots. You can download it from CRAN or there are more features (including formula input and separate colours) in the development version on GitHub: https://github.com/TomKellyGenetics/vioplot

```
devtools::install_github(“TomKellyGenetics/vioplot”)
library(“vioplot”)
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident)

```

To plot several PCs:

```
library(“RColorBrewer”)
colours <- brewer.pal(length(unique(pbmc_small@meta.data$ident)), “Set3”)
par(mfrow = c(1, 2))
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident, main = “PC1”, xlab = “Identity”, col = colours)
vioplot(pbmc_small@dr@pca[,2]~pbmc_small@meta.data$ident , main = “PC2”, xlab = “Identity”, col = colours)
par(mfrow = c(1, 1))

```

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


# Question
Title: How to run proteowizzard on linux using wine?
I am trying to run proteowizzard using wine to be able to convert raw mass spectrometry files to myML format. Following, these instructions\` which seem terribly outdated. I am getting many issues. I am running Ubuntu 16.04.

We want to convert output from a ThermoFisher instrument to mzML or mzXML format.

```
> ~/software/proteowizzard>wine msconvert.exe
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:heap:HeapSetInformation (nil) 1 (nil) 0
> fixme:process:SetProcessShutdownParameters (00000380, 00000000):
> partial stub. err:module:import_dll Library VCOMP110.DLL (which is
> needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\baf2sql_c.dll") not
> found err:module:import_dll Library baf2sql_c.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library VCRUNTIME140.dll (which is needed
> by L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-convert-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-environment-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library
> api-ms-win-crt-filesystem-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library MSVCP140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library VCRUNTIME140.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-runtime-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-string-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-stdio-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-filesystem-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-time-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library
> api-ms-win-crt-convert-l1-1-0.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-environment-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-locale-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-heap-l1-1-0.dll (which is
> needed by L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll")
> not found err:module:import_dll Library api-ms-win-crt-math-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll (which
> is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msparser.dll") not found
> err:module:import_dll Library msparser.dll (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:import_dll Library api-ms-win-crt-utility-l1-1-0.dll
> (which is needed by
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe") not
> found err:module:LdrInitializeThunk Main exe initialization for
> L"Z:\\home\\swacker\\software\\proteowizzard\\msconvert.exe" failed,
> status c0000135

```

# Answer
> 3 votes
Try this Docker container.

All the vendor formats work except T2D, YEP, and FID.

# Answer
> 1 votes
I sympathize ... If it is a Mac just run 'PlayonMac' its a brilliant frontend for Wine without all the bugs and hassles it causes. Its free and available here

---
Tags: mass-spectrometry
---


# Question
Title: Kruskal-Wallis rank sum test output confusion
I am running a multiple comparison using the non-parametric Kruskal Wallis test (K-W), using the ggpubr library and I am a bit confused about the results.

When i just run the KW-test using "base R" the result is different and I am not sure if there is an issue with the data or I am doing something wrong.

In case of ggpubr im getting significant result, i mean the p value is much lower than the critical threshold, but in case of base R this does not occer. I am posting the code and results figure.

Any suggestion would be highly appreciated why results are varying .. Data

```
library(ggpubr)
df1 = read.csv("NEW_RBP/CORPLOT/RBP_DISEASE_CLUSTER/RBP_DISEASE_C1_C4_AVG.txt",header = TRUE,sep = "\t")
head(df1)
my_comparisons <- list( c("HSC", "LSC"), c("LSC", "Blast"), c("HSC", "Blast") )

ex <- melt(df1, id.vars=c("gene"))

head(ex)

ggboxplot(ex, x = "variable", y = "value",ggtheme = theme_bw(base_size = 30),
          color = "variable", palette = "jco", add = "jitter")+
  stat_compare_means(comparisons = my_comparisons)+ # Add pairwise comparisons p-value
  stat_compare_means(label.y = 50)     # Add global p-value

kruskal.test(HSC ~ LSC, data = df1)

    Kruskal-Wallis rank sum test

data:  HSC by LSC
Kruskal-Wallis chi-squared = 416, df = 416, p-value = 0.4908

```

# Answer
The OP says KW is giving different results between "base R" and ggpubr. Background information: KW test is a non-parametric (assumption free) version of 1-way ANOVA.

The analysis you have presented provides a clear result: HSC and Blast are of an approximately equivalent distribution, however LSC is an outlier to both HSC and Blast.

The global value is very highly significant (because of LSC). In other words the global value (HSC, Blast and LSC used in a single analysis) is a much lower probability than the pairwise tests because it will amplify the discrepancy between LSC-HSC and LSC-Blast.

I personally don't see any discrepancy in behaviour within this analysis.

I suspect that HSC and Blast might be significantly different in different analyses because the result is quite low p=0.14, although it isn't borderline (e.g. if p=0.5 I would be happy). The basic issue in my opinion is that non-parametric tests are not particularly powerful and p=0.14 isn't hugely above p=0.05.

If this concerns you, I would examine the HSC and Blast groups using 1-way ANOVA and see if this sheds light on the problem, i.e. the global test will also be rejected by 1-way ANOVA.

1-way ANOVA requires the assumption of normality, but would be supportive of the KW-result. If 1-way ANOVA resulted in HSC and Blast being not significantly different from each other, that is your answer because 1-way ANOVA is a more powerful test (i.e, statistical power), thus it is more likely to reject a relationship between HSC and Blast.
> 0 votes

# Answer
Your problem comes from using an incorrect test, the formula parameter says:

> a formula of the form response ~ group where response gives the data values and group a vector or factor of the corresponding groups.

When you apply the kruskal test, you don't have groups in the LSC variables! You aren't comparing the mean of several groups based on LSC! When you use `stat_compare_means` it is doing a wilcox.test (it hints to it in the help page "a list of additional arguments used for the test method. For example one might use method.args = list(alternative = "greater") for wilcoxon test.").

When you use the wilcox.test we get the same p-value:

```
wilcox.test(df1$HSC, df1$LSC)

    Wilcoxon rank sum test with continuity correction

data:  df1$HSC and df1$LSC
W = 110980, p-value = 4.849e-12
alternative hypothesis: true location shift is not equal to 0

```

In addition, probably you should read about the paired argument, if this two values where measured on the same samples you should use it (but I can't tell it for sure because I don't know what does it mean each variable name).
> 2 votes

# Answer
Just to mention the KW test is appropriate here and not an error as suggested.

Each variable is a compartment and you have 3 compartments. There is no satisfactory equivalent to a non-parametric 2-way ANOVA (possible Friedman test), which is what is being implied.

Wilcoxon will only compare 2 samples, there is no restriction in the number under KW.

If you are comparing alot of samples, KW is much preferred because of type 2 error which will result from Wilcoxon.

I would present both results,KW and Wilcoxon. The only error was that we didn't know what test was run in base R.

Furthermore if you can split each sample into say two compartments, then you would run KW on each compartment (ie run it twice), at this point Wilcoxon becomes erroneous because you would need to run it 6 times at 0.05 critical probability and the chance of a false positive (two data sets are different when they are not) becomes statistically significant.
> 2 votes

---
Tags: r
---


# Question
Title: How to run same command on multiple files?
Hi im doing Variant calling on fastq files for which i have 4000 fq files and the variant calling are done in different 9 steps. Each step generate different files that are output to next steps.

So can anybody help me with this to do it in a loop it pick up names from previous steps for all files. Step 1 command is here:

```
bwa mem -M -R '@RG\tID:ERR040140\tLB:2517707\tPL:ILLUMINA\tPM:HiSeq2000\tSM:ERR040140' ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta ERR040140_* > ERR040140.sam

```

So i made an excel file to break the command in columns so that i can change the file names keeping remaining command same just like this:

So changes are occurring 5th, 7th and 9th column. which are:

```
tID:ERR040140
tLB:2517707
tSM:ERR040140
ERR040140_*
ERR040140.sam

```

Can anyone help me with loop or something?

Thank you!

# Answer
Another solution for repeating one command with different parameters is to use `parallel`.

Create a tab-delimited file `samples.txt`, which contains the sample names and the information for your read groups, e.g.:

```
Sample1 Lib1
Sample2 Lib2
Sample3 Lib3

```

Then you can use `parallel` like this:

```
parallel --colsep "\t" 'bwa mem -M -R "@RG\tID:{1}\tLB:{2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:{1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta {1}_* > {1}.sam' :::: samples.txt

```

For each line in `samples.txt` a new process is started, replacing `{1}` and `{2}` with the values in column 1 and 2 of this line.

fin swimmer
> 10 votes

# Answer
Something like this might do it:

```
#!/bin/bash

# make an array to hold your list of samples/base filenames
samples=( S1 S2 S3 S4 S5 )

# make a loop to go through the sample IDs
for sample in "${samples[@]}"
do
  :
  # do something with each sample/file, e.g. make a directory for it
  mkdir "${sample}"_home

  # do something else with the sample, e.g. unzip the fastq.gz file
  gunzip "${sample}".fq.gz > "${sample}"_home/"${sample}".fq
done

```

One disadvantage with this is that it will process the files one after the other. If you have a cluster, you might prefer to run them all at the same time to make the whole thing quicker. If you have a SLURM cluster, you can use job array to do something similar.

If you have your list of samples in a text file, you could try replacing the samples array command in the script above with something like this:

```
samples=( `cat samples_list.txt` )

```

This will only work if the samples are listed on a single line, space-separated. You may need to process your samples file further to get them in this format, or use the bash readarray command.
> 8 votes

# Answer
if you have a repeating workflow, I strongly recommend to have a look at workflow management systems like snakemake. I've also wrote a little tutorial on biostars about this topic, which might be useful for you.

fin swimmer
> 5 votes

# Answer
You can create a script to pass arguments to as follows:

```
#! /bin/bash

bwa mem -M -R "@RG\tID:${1}\tLB:${2}\tPL:ILLUMINA\tPM:HiSeq2000\tSM:${1}" ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta ${1}_* > ${1}.sam

```

Now if you run `bash script.sh sample lib` then `sample` and `lib` will be the script as arguments `${1}` and `${2}` respectively. Then all you need is a loop to pass arguments to the script:

```
for sample in sample1 sample2 sample3
   do
   for lib in lib1 lib2 lib3
     do
     bash script.sh $sample $lib
     done
   done

```

You can use `nohup` to run each sample in the background and proceed with next one. This will in essence, submit each command as a separate job to run on a different core. These will print to separate log files sample.out instead of to the terminal. It will also continue to run, even if you disconnect from SSH or close the terminal.

```
for sample in sample1 sample2 sample3
   do
   for lib in lib1 lib2 lib3
     do
     nohup bash script.sh $sample $lib > ${sample}.out &
     done
   done

```

You should only do this for single-threaded processes. For multi-threaded processes or memory-intensive jobs, this will demand a lot from the machine. Even if this isn’t a shared server, it will perform poorly if you run out of memory and have to “swap” memory by writing to disk. So it’s best to submit them one-by-one (above) or use a cluster with a job scheduler.
> 0 votes

---
Tags: variant-calling, linux
---


# Question
Title: Can blat use more than one core/CPU to speed up the alignment?
I am using BLAT to align two versions of the genome of *C. elegans*. I can see in the Activity Monitor of my Mac Book Pro High Sierra that `blat` is using 100% of a CPU. However, is this programme able to use more than one core / CPU to speed up the alignment?

I cannot find such option in the help page of the tool here.

I installed blat using `conda install -c bioconda blat`.

I am using the latest version of `blat`, i.e. 36.

# Answer
> 6 votes
BLAT can only use one CPU. It is actually not the right tool for full-genome alignment. For "two versions" of the same species, MUMmer and minimap2 are orders of magnitude faster and probably give better alignment.

EDIT (moving comment to answer): OP comments that the purpose of this alignment is for lifting over annotations using the UCSC PSL-based pipeline. So far there are no converters for minimap2's SAM/PAF or MUMmer's delta formats. With the UCSC pipeline, you have to write a converter by yourself. An alternative approach is to use minimap2's liftover. The minimap2 cookbook gives an example:

```
minimap2 -cx asm5 --cs ecoli_ref.fa ecoli_canu.fa > ecoli_canu.paf
echo -e 'tig00000001\t200000\t300000' | paftools.js liftover ecoli_canu.paf -

```

Note that this pipeline is not as mature as UCSC liftOver and is not as scalable.

---
Tags: sequence-alignment, software-usage, blat
---


# Question
Title: Experimental Design for Differential expreression analysis
I have a Normal esophageal Fibroblasts (NOFs) cultured in DMEM media; The same NOF also have been cultured with a tumor sample from a patient named 005 on DMEM media; I have also Cancer Associated Fibroblasts (CAFs) cultured on DMEM, on DMEM + tumor from a patient named as 005 and DMEM + tumor from a patient named 036 something like this picture;

I want to know if NOF and CAF are different in terms of gene expression but I am not sure about design for DESeq2 or Edger and ANOVA.

By @Devon's help I have done so

```
> mycols
   cell_type co_cultured
G2       NOF          NO
G3       NOF         YES
G4       CAF          NO
G5       CAF         YES
G6       CAF         YES

dds=DESeqDataSetFromMatrix(countData = NOFCAF,colData = mycols, design        =~ cell_type*co_cultured)

```

Now, I am wondering what would be a suitable contrast for results function in DESeq2

```
> resultsNames(dds)
[1] "Intercept"                   "cell_type_NOF_vs_CAF"        "co_cultured_YES_vs_NO"
[4] "cell_typeNOF.co_culturedYES"

> res=results(dds)
> res
log2 fold change (MLE): cell typeNOF.co culturedYES
Wald test p-value: cell typeNOF.co culturedYES
DataFrame with 2545 rows and 6 columns

```

Thanks a lot for correcting me

# Answer
> 3 votes
You're really running out of degrees of freedom and the only actual replicate is a single sample from another patient, so please take the results with a huge grain of salt.

You can use a design of `~cellType * coculture` where `cellType` has levels `NOF` and `CAF` and `coculture` has levels `yes` and `no` (ideally you'd block by patient, but that's not an option for you).

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


# Question
Title: Can structural - equation modeling be used in Bioinformatics?
I am a bioinformatics student and I was learning the basics of structural equation modeling (SEM).

I was wondering, are there any applications of SEM on bioinformatics? Are there any strong sides for use of structural equation models using data on bioinformatics?

I did a quick search and found almost nothing. While it seemed like a learning algorithm, I can't say I figured it out completely.

# Answer
> 1 votes
I have came across them in gene regulatory network analysis before. see here and here. However I haven't seen them in other contexts, but that's not to say they don't exist.

# Answer
> 1 votes
SEM specialises in identifying 'latent variables'.

SEM are extensively used for questionnaire data. In questionnaire data, latent variables are the 'psychological' factors behind given (human) behaviours.

Applying SEM to explicit molecular mechanisms doesn't make much sense unless there is a hidden component in a mechanism that has not been identified. However, even here the SEM will not tell you what the 'latent variable' represents.

Thus I don't believe there is a good application to bioinformatics, outside human phenotypes derived from questionnaire data, which might be mapped on biological variables, but I could be wrong.

---
Tags: statistics, algorithms, structural-variation
---


# Question
Title: sva for RNA-Seq data without known phenotype
I have been working on RNA-Seq data from two different cohorts, and they show very strong batch effect (~35% variance explained by 1st component in PCA). Since I am trying to do a class discovery from a data set with the subtype of only some samples are known, the only methods I have been using are ComBat and pSVA from SVA package. I would like to compare them with more methods to make the clustering result more trust-worthy.

I found that some papers used sva package, but not citing the paper on pSVA, meaning perhaps there are some other ways to use sva for this type of data that I am not aware of. Is there a method to construct a null model and a full model for sva only by knowing the batch the samples are from?

Any other suggestion of batch effect correction tools for a data set with no known phenotype/subtype are also welcome.

# Answer
> 1 votes
You can just remove the first component from the dataset by setting the first eigenvalue to 0 in the diagonal matrix and then multiplying the SVD matrices.

I am not sure which R code you are using to estimate the components (I have a comparison of PCA code in R on https://github.com/aedin/ODSC\_2018/tree/master/PCA\_vignette)

In ade4 R package, there is a function `reconst()` which reconstructs data from a selection of components: https://rdrr.io/rforge/ade4/man/reconst.html

Happy to provide code to do this in other PCA versions if needed Aedin Culhane

---
Tags: rna-seq, clustering, batch-effects, combat, sva
---


title: GDC data portal cancer subtype

text:
# Question
Title: GDC data portal cancer subtype
I need to use Gene Expression quantification data (using HTSeq-FPKM-UQ) for breast cancer data on the gdc portal (https://portal.gdc.cancer.gov/repository) to classify breast cancer subtypes (luminal A, luminal B, HER2-like, basal-like).

However, since I don't know the portal, I'm having problems in associate the EXP data and the label.

I could get the label either by using the old TCGA-barcode and use a dictionary I already have from an article (but I don't know how to get the barcode from the portal / API) or by using clinical data (but I don't know how to find metadata that link the EXP and the clinical data).

Is there anyone more used to the portal who could help me?

# Answer
> 1 votes
I personally find the GDC portal not very user-friendly. However, some of the data is available from other sources. For example, you can download TCGA-BRCA FPKMs and phenotypes from UCSC Xena here: https://xenabrowser.net/datapages/?cohort=GDC%20TCGA%20Breast%20Cancer%20(BRCA)

All the sample info is already combined and the labels are the classic TCGA barcodes.

# Answer
> 1 votes
These data are in Experiment Hub in Bioconductor

```
library(curatedTCGAData)
brcaRNA <- curatedTCGAData("BRCA", "RNASeq2GeneNorm", dry.run = FALSE)

```

---
Tags: public-databases, cancer, data-download, data-retrieval
---


# Question
Title: Using multidimensional scaling to visualize protein sequences by functionality
I have a multiple sequence alignment (MSA) of protein sequences, with which I have performed multidimensional scaling to visualize their clustering. Using R, I used the `seqinr` package to create an identity matrix (using `dist.alignment`). Is there a similar method that is more "generalizable"? Specifically, I'm interested in classifying the amino acids by some functional characteristics, such as polarity, instead of having the alignment list the amino acids by their identity (A, G, I, etc.).

Here's a very short set of fake input data, with columns 2-4 representing amino acid residue positions:

```
a <- c("Strain 1", "Strain 2", "Strain 3", "Strain 4")
b <- c("Polar_hydrox", "Nonpolar", "Nonpolar", "Polar_hydrox")
c <- c("Nonpolar", "Nonpolar", "Nonpolar", "Nonpolar")
d <- c("Nonpolar_cyclic", "Nonpolar_cyclic", "Polar_amide", "Polar_amide")
e <- c("Polar_amide", "Acidic", "Polar_amide", "Acidic")
x <- data.frame(a, b, c, d, e)

```

Is what I'm asking a feasible/correct use of MDS? Or should I be looking into a different method entirely? I attempted to use Gower's formula to create a dissimilarity matrix on the original dataset (i.e. using the actual amino acid identities instead of the functional groupings), but compared to what I got using `dist.alignment`, the MDS plots look different.

Edit: Gower's formula and the `dist.alignment` function are obviously calculating the distance matrices differently, but I don't know the details on the different methods. Would using Gower's formula be a good approach for calculating distance matrices for the example dataset?

Edit: clarified functions used for calculating distance matrix.

# Answer
Some of the first MSA analysis, eg of codon usage of bacteria were performed using matrix factorization (see old papers from Des Higgins). His group also published discriminative methods for finding functional residues using supervised correspondence analysis

> Wallace IM1, Higgins DG BMC Bioinformatics. 2007 Apr 23;8:135. Supervised multivariate analysis of sequence groups to identify specificity determining residues.

The method is available in the R package ade4 and Bioconductor package made4

Also have a look at the phylogenetic R package adephylo and several related packages
> 2 votes

# Answer
I personally don't think this will work well for some variables, primarily because you are not mapping conformational epitopes, where amino acids will cluster in physical space by their charge. In particular, the charge is critical to retaining this conformational structure.

However, hydrophobic/hydrophilic characters are far more amenable to linear style epitopes that this analysis is primed towards.

I don't recommend a multi-dimensional array for this analysis, because you have loads of information, but rather a 'dictionary' (Python), hash (Perl) (or 'object' \[any\]) or XML, which retains all the information for a given sequence and then strips this in preparation for your multivariate statistic.
> 1 votes

---
Tags: r, clustering
---


# Question
Title: How BLOSUM Matrix is constructed and calculated
I would like to ask how BLOSUM matrix is constructed and calculated ? I read the wikipedia and I still do not understand it. I do not understand the mathematical calculations as I have low knowledge about logarithms.

Maybe could someone give a simple example and then a real-scenario example of how the BLOSUM matrix can be used and calculated, maybe in relation to the Blast algorithm where it uses the BLOSUM scoring matrix to determine high scoring words for each word in the query sequence ?

Also, why are they different types of BLOSUM matrices like BLOSUM62 and so on ?

# Answer
> 3 votes
Sean Eddy wrote an excellent explanation of the BLOSUM Matrix in his Nature Biotechnology paper Where did the BLOSUM62 alignment score matrix come from? He also gives real examples of how one would calculate the score for different amino acid pairs and where the numbers come from. Further references given at the end of the paper might also be useful reading, specifically: Amino acid substitution matrices from an information theoretic perspective, Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes, and Amino acid substitution matrices from protein blocks.

# Answer
> 2 votes
BLOSUM stands for BLOcks SUbstitution Matrices (Henikoff & Henikoff, 1992), and were created by observing substitution frequencies in local ungapped multiple sequence alignments. The score reflects the chance (log-odds) one amino acid is substituted for another in a set of protein multiple sequence alignments. BLOSUM scoring matrices are normally followed by a number eg BLOSUM62. The number is % sequence identity between the sequences in the multiple sequence alignment (MSA) used to create the score matrix. eg BLOSUM62 matrices were created from multiple sequence alignments with blocks that shared 62% identity. Note because these are ungapped alignments, the multiple sequence alignments were generated from conserved protein families. With the older PAM matrices, the number is opposite, a higher number is greater distance. PAM matrices reflect evolutionarily distant in global (aka gapped) alignments of the full-length sequences that include both conserved and diverged regions. See Wikipedia https://en.wikipedia.org/wiki/BLOSUM. There are several newer substitution matrices. The R package SubVis allows you to compare/visualize different substitution matrices.

# Answer
> 1 votes
Amino acid matrices are contructed by maximum likelihood and easily performed using PAML. You present a known tree and an amino acid alignment and ask what amino acid matrix would maximise the likelihood. Hence all amino acid matrices using this old school approach are dependent on the alignment composition. JTT was generally regarded as obtaining the highest likelihood. However, I reitereate that things have moved on a lot since this 2-D reversible amino acid substitution matrix was devised.

---
Tags: sequence-alignment, blast
---


title: Is there any R package for biodbnet?

text:
# Question
Title: Is there any R package for biodbnet?
I want to convert Ensembl ids to Gene\_symbols. I used biomaRt package but all Ensembl ids are not seen in the results.

Examples of problematic IDs are: ENSG00000255414, ENSG00000274358, ENSG00000235565, ENSG00000278309. I gave them in bioDBnet, which seems to handle them.

I would like to know whether there is any R package for bioDBnet or not? If so, could you please tell how to use it?

# Answer
I couldn't find a package neither in Github nor in Bioconductor. Neither could I find any paper about a package linking to bioDBnet. So I don't believe there is any (R) package at the moment which links to biodbnet in order to retrieve the IDS.

However it has a couple of APIs so you could build one or use them to retrieve the information you want.
> 4 votes

# Answer
```
library(EnsDb.Hsapiens.v86)
library(ensembldb)
genes(EnsDb.Hsapiens.v86)
genes(EnsDb.Hsapiens.v86, filter=GeneIdFilter('ENSG00000235857'))
transcripts(EnsDb.Hsapiens.v86, filter = GeneIdFilter("ENSG00000184895"),
              order.by = "tx_seq_start")

```

(The ids you listed are older than v86). See the vignette for the ensembldb Bioconductor package.
> 0 votes

---
Tags: r, bioconductor, ensembl
---