cannin/biostars_qa huggingface dataset:
Question: I have a MeatG library, which was sequenced by two times. The first time, we got unpaired reads (i.e., unpaired.fastq); while we got interlaced forward and reverse paired-end reads (i.e., interlaced-paired.fastq) for the 2nd time. Since they are from the same library, I assembled them together and got the assembled scaffolds. Now I would like to align both of them (unpaired and interlaced paired reads) simultaneously to the assembled contigs using Bowtie2 version 2.3.4.1 to get one .sam file. Can Bowtie2 do this? Any suggestions are appreciated!

I tried the following two scripts, the results for both were same. So for the first script, only the paired reads are mapped.

bowtie2 -q --phred33 --end-to-end --sensitive -p 12 -I 0 -X 2000 --no-unal -x bowtie2-db/CBIW --interleaved interlaced-paired.fastq -U unpaired.fastq | samtools view -Sb - > output.bam

bowtie2 -q --phred33 --end-to-end --sensitive -p 12 -I 0 -X 2000 --no-unal -x bowtie2-db/CBIW --interleaved interlaced-paired.fastq | samtools view -Sb - > output.bam

**See the align result:**

Building a SMALL index

35388893 reads; of these:

35388893 (100.00%) were paired; of these:

11580370 (32.72%) aligned concordantly 0 times

23621928 (66.75%) aligned concordantly exactly 1 time

186595 (0.53%) aligned concordantly >1 times

11580370 pairs aligned concordantly 0 times; of these:

101972 (0.88%) aligned discordantly 1 time

11478398 pairs aligned 0 times concordantly or discordantly; of these:

22956796 mates make up the pairs; of these:

22513992 (98.07%) aligned 0 times

432094 (1.88%) aligned exactly 1 time

10710 (0.05%) aligned >1 times

68.19% overall alignment rate

Answer: I don't think Bowtie2 can map paired- and single-end reads simultaneously, but you can merge both mappings with samtools. BWA can map map a mix of interleaved paired reads and single reads, if you concatenate both files and keep single reads at the end.

Use `samtools sort` to sort each bam by position, then `samtools merge` to merge the sorted bams into one.

If having a sorted bam is not important, you can just use `samtools cat`.
Question: I have files with the 2504 individuals of the 1000 genomes project, and I want to filter by population. I did the following for the first population (ACB):

plink --file all1000gen --keep indACB.txt --make-bed --out all1000genACB
but it gives back the following error:

Error: Line 1 of --keep file has fewer tokens than expected.
my indACB.txt file looks like this:

head indACB.txt
HG01879
HG01880
HG01882
HG01883
HG01885
HG01886
HG01889
HG01890
HG01894
HG01896
which I made (por each population, using grep) from the population information file that's available in the 1000 genomes page, which has a two times the individual ID (first two columns) and one with the population name, as shown:

head indpop2.txt
HG00096 HG00096 GBR
HG00097 HG00097 GBR
HG00099 HG00099 GBR
HG00100 HG00100 GBR
HG00101 HG00101 GBR
HG00102 HG00102 GBR
HG00103 HG00103 GBR
HG00105 HG00105 GBR
HG00106 HG00106 GBR
HG00107 HG00107 GBR
I think there's a problem with my --keep file, but I'm not sure what's the wanted structure of the txt file.

I also tried greping ACB individuals from indpop2.txt , so the new indACB.txt file looks like this:

head indACB2.txt
HG01879 HG01879 ACB
HG01880 HG01880 ACB
HG01882 HG01882 ACB
HG01883 HG01883 ACB
HG01885 HG01885 ACB
HG01886 HG01886 ACB
HG01889 HG01889 ACB
HG01890 HG01890 ACB
HG01894 HG01894 ACB
HG01896 HG01896 ACB
But it yields the following error:

plink --file allconcat39 --keep indACB2.txt --make-bed --out allconcat43ACB

Error: No people remaining after --keep.

Answer: from https://www.cog-genomics.org/plink/1.9/filter

> --keep accepts a space/tab-delimited text file with family IDs in the first column and within-family IDs in the second column, and removes
> all unlisted samples from the current analysis.

So like this:

HG00096 HG00096
HG00097 HG00097
HG00099 HG00099
HG00100 HG00100
HG00101 HG00101
HG00102 HG00102
HG00103 HG00103
HG00105 HG00105
HG00106 HG00106
HG00107 HG00107
Question: I have two fasta files, with the same headers/names for the sequences but different sequences.
I would like to combine them into one file, so that each sequence has the same name but is a combination of both sequences.
My preferred language is bash script, but I'm open to other suggestions.
thanks.

Answer: A solution using [seqkit](http://bioinf.shenwei.me/seqkit/), [csvtk](http://bioinf.shenwei.me/csvtk/) and shell `sed`.

Sample files (**not in same order**, **can be multiple lines**):

$ cat 1.fa
>seq1
aaa
aa
>seq2
ccc
cc
>seq3
ggg
gg

$ cat 2.fa
>seq3
TTT
TT
>seq2
GGG
GG
>seq1
CCC
CC

Just one command:

$ seqkit concat 1.fa 2.fa
>seq1
aaaaaCCCCC
>seq2
cccccGGGGG
>seq3
gggggTTTTT

**Step 1**. Convert FASTA to tab-delimited (3 columns, the 3rd column is blank (no quality for FASTA)) file:

$ seqkit fx2tab 1.fa > 1.fa.tsv
$ seqkit fx2tab 2.fa > 2.fa.tsv

$ cat -A 1.fa.tsv
seq1^Iaaaaa^I$
seq2^Iccccc^I$
seq3^Iggggg^I$

**Step 2**. Merge two table files:

$ csvtk join -H -t 1.fa.tsv 2.fa.tsv | cat -A
seq1^Iaaaaa^I^ICCCCC^I$
seq2^Iccccc^I^IGGGGG^I$
seq3^Iggggg^I^ITTTTT^I$

**Step 3**. Note that there are two TAB between the two sequences, so we can remove them to join the sequences

$ csvtk join -H -t 1.fa.tsv 2.fa.tsv | sed 's/\t\t//'
seq1 aaaaaCCCCC
seq2 cccccGGGGG
seq3 gggggTTTTT

**Step 4**. Convert tab-delimited file back to FASTA file:

$ csvtk join -H -t 1.fa.tsv 2.fa.tsv | sed 's/\t\t//' | seqkit tab2fx
>seq1
aaaaaCCCCC
>seq2
cccccGGGGG
>seq3
gggggTTTTT

**All in one command**:

$ csvtk join -H -t <(seqkit fx2tab 1.fa) <(seqkit fx2tab 2.fa) | sed 's/\t\t//' | seqkit tab2fx
Question: Hello everyone,

I need some help/ideas in calculating average depth per SNP across samples. I have already generated depth files for each sample (around 18000 samples in total) using `samtools depth input.bam > output.txt` and I have merged all the `output.txt's` into 1 txt.gz. The input file is sorted by position and it looks like this.

sample1 17 123456 1000
sample2 17 123456 500
sample3 17 123456 10
sampleN 17 123456 50
sample1 17 123457 1100
sampleN 17 123457 50

so first column contains sample names second column is chromosome number third column is the snp position and the last column is depth. I am just wondering if there is quick way to get average depth per SNP across samples as my input file is huge (around 38gb)? or if there is an alternate way?

Looking forward to your suggestions and ideas.

Regards,

Answer: You can try with datamash:

$ datamash -g 2,3 mean 4
Question: Hi,

I want to perform GSEA analysis with `clusterProfiler` package and I am following [this post][1] . But I have few doubts (and be patient, I am a biologist):

1. Is it better to perform GSEA using the Ensembl ID or the gene symbol?

I did it with both Ensembl & gene names just to see if I was getting the same results but surprisingly I got slightly different results. I am not sure which one should I choose since I don't want to be biased based on what would fit/support better my hypothesis (which is the one run on the gene symbol).

2. in the first figure (dotoplot) of the website I posted, the enriched terms are divided by activated and suppressed which reflect whether the enrichment score (ES) is positive or negative. A positive ES indicates gene set enrichment at the top of the ranked list whereas a negative ES indicates gene set enrichment at the bottom of the ranked list. but what does it mean practically? What does it tell me that e.g. cytosol ribosome (in the figure) has negative ES so it is at the bottom of the ranked list?

thank you in advanced

Camilla
[1]: https://learn.gencore.bio.nyu.edu/rna-seq-analysis/gene-set-enrichment-analysis/

Answer: Hi Camilla,

I just recently started working with this package but i will do my best to answer your question to the best of my current knowledge.
About your first question:
> Is it better to perform GSEA using the Ensembl ID or the gene symbol?

I personally would stick to ENSEMBL IDs as it is very often the case that multiple ENSEMBL IDs have the same gene symbol. I think this might also explain why your results would differ using geneSymbols vs ENSEMBL IDs. You can simply check how that on your input data using R: `length(unique(YOURDATA$ENSEMBLID))`& `length(unique(YOURDATA$SYMBOL))`. This will give you an idea on how much your input varies in size for the GSEA analysis. Ultimately when there are fewer genes to compare with in the enrichment model this will influence the outcome. So I would say using the dataset with the maximum of information (I assume that would be the ENSEMBL IDs) would be the less biased approach.

> [...] the ranked list. [...] what does it mean practically?

As far as I understand the concept of the running enrichment score and the approach of clusterProfiler you have to provide a gene list where the effect size (i.e. log2(fc)) is sorted from high to low. Hence, your input data is already ranked in a way that up regulated genes (positive log2(fc)) have lower ranks than your downregulated genes. This explains why negative ES have a lower ranks compared to those with a positive ES.

I apologize if my current explanation is a little bit cryptic but you can run the `gseaplot()` function, which displays you the running ES over the position in the ranked list of genes for a particular pathway in your results list. Where this plot reaches it's maximum you see a dashed red line. I think where this line intersects with the y-axis it determines the rank position of this particular pathway. Please [check here][1] for details.

I hope this helps :)

[1]: http://yulab-smu.top/clusterProfiler-book/chapter12.html#running-score-and-preranked-list-of-gsea-result
Question: I'm running BLAT with this command:

blat genome.fasta query.fasta blat_90.csv -minIdentity=90 -noHead

And I'm surprised that I've counted the results and order by count.

import pandas as pd
import sys
df_blat = pd.read_csv('blat_90.csv', index_col=False, sep='\t')
s = df_blat.qName.value_counts()
count = 1
for k in s.head(100).index.tolist():
count += 1
print(s[k],k)

The interesting thing is that the max results per query is 672 and this number is repeated about 30 times. Then I have queries that were repeated 671 times again many times.

672 MITE_T_71978|chr3D|340455280|340455443|AT|175|F2975
672 MITE_T_97305|chr7B|194283260|194283371|AT|29|F4298
672 MITE_T_110543|chr5D|34445518|34445608|TA|94|F5023
672 MITE_T_72475|chr3D|503092467|503092630|AT|174|F2988
...

So this makes me wonder if there's a limit in the number of times a query is search. I've repeated this script with no different results.

UPDATE:
I've notice that running only one query with each wheat chromosome in a separate thread, the max results per search is 32.

This is how the results look like when I search only one sequence agains each chromosome at the time:

wc -l b*
32 b1A.csv
32 b1B.csv
32 b1D.csv
32 b2A.csv
32 b2B.csv
...

Max 32 hits per search (accordingly 32 x 21 chromosomes is 672)

There must be a limit in the code but I'm still not able to find it

Answer: As responded by UCSC itself:

> The 32 results are indeed a consistent limitation where BLAT only
> allows 16 results per chromosome and per strand. BLAT is unable to
> exhaustively return results for queries with repetitive sequences that
> map to many places in the genome. There may be many other sequence
> aligning programs available that do not have such limits, but we do
> not have any specific recommendations about which to use.
>
Question: I have sought previous posts (post1, post2) on how to call somatic variants and it seems the general practice is to intersect multiple callers to insure that low VAF mutations are being called.

My approach was to samtools, MuTect2, SomaticSniper, VarScan2 but I found an interesting post saying that as long as read placements are perfect, any caller suffices (even samtool mpileup). I should mention that I am working with RNA seq of cancer samples (matched with normal).

In general, my view is as long as read placement is perfect, even naive methods work sufficiently well... To me, the simplest yet most effective strategy is to use two distinct alignment algorithms, such as bwa and bwa-sw, which have distinct error modes. You only consider mutations shared between the two alignments... Another complication is structural variations, in which I am less experienced. In some sense, false mutations caused by structural variations are still indication of something different between normal and tumor... In all, I think you do not need to worry about which software to use for detecting somatic mutations - anything reasonable is fine. You should pay more attention to mismapping and structural variations.

First, does read placement mean how well aligners align our samples to the reference? Does working with RNA-Seq introduce higher error rates in read placements? Has the consensus changed or is intersecting multiple callers still recommended? Thank you very much for your help.

Answer: If someone tells you that somatic mutation calling is easy or a solved problem, they have never really tried to do somatic mutation calling.

There are a host of issue to contend with - sequencing artifacts, problems with the reference, differential coverage, (and yes, mismappings are common!), etc. Your approach of using mutliple callers seems sensible. The tricky part is figuring out how to combine them. Straight intersections will give you high specificity, but low sensitivity. Unioning the three will result in the opposite. A more nuanced approach has been been explored by recent tools like [somaticseq](http://www.genomebiology.com/2015/16/1/197). I haven't used that one in particular, but I am convinced that an approach of that nature is most likely to succeed.
Question: Hi,

Am looking for a fast way to parse through 2 large lists with millions of elements and identify the ones shared by both the lists, exclusive to list1, exclusive to list2. For smaller lists, perl hashes or the venny tool are useful, but I needed to do the same with these huge lists.

FWIW, these lists are actually read-ID from NGS data. The reads were aligned using 2 different approaches, and I wish to investigate of the ~ 30 million reads mapped by both approaches, how many are common, how many and which ones are aligned using only one approach and not the other. I mention this in case there is something in bedtools, bamtools or the like that might be relevant here.

thanks!

EDIT: Thanks for all the responses. I got bogged by other things and never got the chance to check on these. Will work my way through and add further comments/notes..

Answer: Assuming you are referring to IDs and you have lists that are structured something like the following:

$ cat list-1.txt
foo
bar
baz
zab
oof

$ cat list-2.txt
zab
baz
zit

You can first sort the two files:

$ sort list-1.txt > list-1.sorted.txt
$ sort list-2.txt > list-2.sorted.txt

$ cat list-1.sorted.txt
bar
baz
foo
oof
zab

$ cat list-2.sorted.txt
baz
zab
zit

Now, you can use join to find the elements that are common to both sets:

$ join -1 1 -2 1 list-1.sorted.txt list-2.sorted.txt
baz
zab

Use the -v option to find those elements that are exclusive to set 1:

join -v 1 -1 1 -2 1 list-1.sorted.txt list-2.sorted.txt
bar
foo
oof

Exclusive to set 2:

join -v 2 -1 1 -2 1 list-1.sorted.txt list-2.sorted.txt
zit

Or, you can use the comm command to do it all in one step. The first column are elements exclusive to set 1, the second column is elements that are exclusive to set 2, and the third column is elements that are common to the two sets. join is more flexible, but in simple cases like this, comm can be quite useful.

$ comm list-1.sorted.txt list-2.sorted.txt
bar
baz
foo
oof
zab
zit
Question: I have genotyping data coming from GenomeStudio (Illumina) software from a SNP-chip, and I would like to convert this data to VCF. Does anyone know of a tool to do this? or should I write my own script to do this?

Here is how the GenomeStudio output file is formatted:

Index Name Address Chr Position GenTrain Score Sample1 Sample2
1 cnvi0000001 124680304 5 164388439 0 AA AB
...

Answer: Genomestudio has a [PLINK export plugin][1], and you can use [PLINK/SEQ][2] to convert that format to a VCF.

[1]: http://support.illumina.com/downloads/plink_input_report_plug-in_v213.ilmn
[2]: https://www.researchgate.net/go.Deref.html?url=http%3A%2F%2Fatgu.mgh.harvard.edu%2Fplinkseq%2F
Question: I'm a Computer Science undergraduate writing my thesis which is about analyzing open source bioinformatics projects and extracting standalone modules. I need those projects to be written in Java and preferably be stable or mature. So far I've found a bunch of projects via SourceForge and I was wondering if there are any notable ones that are not hosted there. A different repository that has a bioinformatics category would be helpful too.

Here's a list of the most promising projects I have so far:

- Jenetics
- JGap
- Jmol
- Juicebox
- Picard
- RDP Classifier
- Shim
- Tassel
- The Chemistry Development Kit
- Toxfree
- TreeView
- VarScan

Thanks in advance.

Answer: - [MHAP][1] (part of [Canu][2] long read assembler)
- The server part of [QTLNetMiner][3] is written in Java https://github.com/KeywanHP/QTLNetMiner
- PacBio's SMRT Analysis is partially written in Java, but I don't know whether it fits your requisites, it seems to be open source but I can't find a license
- [JAligner][4]
- [Jalview][5]

[1]: https://github.com/marbl/MHAP
[2]: https://github.com/marbl/canu/tree/master/src/mhap
[3]: https://ondex.rothamsted.ac.uk/QTLNetMiner/
[4]: http://jaligner.sourceforge.net/
[5]: http://www.jalview.org/
Question: I usually use `bedtools intersect` to find overlapping regions of bed files, but it seems like this tool can only output overlap between a pair of files. I need something that can do the following from many bed files and **only** report regions contained in all of the bed files.

Example input from 4 separate bed files:

```
chr1 50 100
chr1 60 120
chr1 30 90
chr1 50 90
```

Desired output:

```
chr1 60 90
```

Any tools for this? Maybe I should just `cat` all the bed files together and merge them?

Answer: To intersect intervals from 1 to `N` files, simply use `bedops --intersect`:

$ bedops --intersect A.bed B.bed ... N.bed > answer.bed

More details from the documentation: http://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#intersect-i-intersect

This requires sorted BED files. You can use [BEDOPS `sort-bed`][1] to do this quickly.

I think the other toolkit suggestions now require sorted BED files, as well.

[1]: http://bedops.readthedocs.io/en/latest/content/reference/file-management/sorting/sort-bed.html
Question: The tool is at: http://bedtools.readthedocs.org/en/latest/content/tools/coverage.html

I want to get base coverage for specific regions (I'll provide a bed file). However, the output is not in a standard BedGraph format that I can export into R for visualisation.

Q: How do I draw density plot of coverage over the regions from the outputs generated by the coverage tool?

Q: Alternatively, how to convert the non-standard outputs to a proper BedGraph? For example, the Sushi R package takes a bedgraph input file.

Answer: You should also consider bedtools genomecov, with the `-bg` option (bg being bedgraph)

http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html

Edit: After reading comments, and after some experimenting, what I think you want to do is:

```
bedtools intersect -a myBam.bam -b myRegions.bed > intersected.bam
bedtools genomecov -trackline -bg -ibam intersected.bam > intersected.bedgraph
```

Bedtools intersect will give only bam reads that overlap your regions (in myRegions.bed)

Then genomecov will draw a bedgraph file (coverage) only for the read depth for your filtered regions.
Question: Is it possible to search for orthologs/inparalogs only for specific sequences in the query genome?

I'm aware that searching for orthologs requires complete proteomes. I would like to avoid a step where OMA calculates all-versus-all sequence comparisons of all proteins from query and subject genomes. In fact, I don't have to BLAST all proteins in a proteome when I'm interested only in one protein. I was wondering whether OMA can be used this way? For example, [Reciprocal Smallest Distance (RSD)][1] approach lets users to find orthologs for specific sequences in the query genome.

[1]: https://github.com/todddeluca/reciprocal_smallest_distance

Answer: If your specific sequence is part of a genome in the OMA database, it is of course possible to retrieve the inferred orthologs for this sequence only. For this, use the search function with an identifier or the sequence itself.

If the sequence is not in OMA, you can export a few closely related genomes and run the OMA standalone on them and your sequence. Assuming that the genomes are somewhat complete, the orthologs inferred for your sequence should be fairly reliable. Thanks to the reciprocal criterion (roughly the RSD criterion you mention) and the use of "witnesses of non-orthology", things only start to fall apart when most genomes are very incomplete.
Question: Hello,

Sorry in advance for my English.

I want extract all positions of non-sequenced nucleotides for each chromosomes from a genome sequence file (Fasta).

To do that, I have a fasta file with the genome sequences, like this:

```
>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNtaaattgttt
taaattgtttctgtttgcagttgacatgatcttatatatagaaaacacca
ataactctgccaaaaaatttagaattcataaatgaatttagtaaagttgc
```

I want find all N positions and obtain this position in bed format, e.g.:

chr1 1 200 ...

I didn't found how to do that in Bedtools.

If you have any idea, could you help me?

Thank you

Answer: Well, this should not be a problematic task if I understood the question correctly. All you need to do is slide along your input fasta file and count the swithches N-> something and something to ->N

perl -ne 'chomp;if( />(.*)/){$head = $1; $i=0; next};@a=split("",$_); foreach(@a){$i++;if($_ eq "N" && $s ==0 ){print "$head\t$i"; $s =1}elsif($s==1 && $_ ne "N"){print "\t$i\n";$s=0}}' infile.fa

This should do the trick.

cheers

mxs
Question: Hello,

I've done a multi-sample SNP calling with samtools mpileup, and my results are in a vcf file.

I would like to keep only SNPs detected on all of my samples (without any sample with GT flag = 0/0 ,(homozygote for the reference allele)).

Does anyone know a tool that could do it?

Thanks!

Answer: I usually choose vcftools or bcftools to deal with vcf files, but in your case a simple grep would do:

grep -v "0[/|]0" <in.vcf >out.vcf
Question: Hi,

How can I remove duplicate rows from the output of my toptableF? I have tried the following without success:

if (unique(fit2$genes)){
Sigenes.F <- topTableF(fit2, number=100, genelist=fit2$genes, adjust.method="bonferroni",
sort.by="F")
}

Thanks

Answer: Read ?duplicated:

> x = c(1:10,1:10)
> x
[1] 1 2 3 4 5 6 7 8 9 10 1 2 3 4 5 6 7 8 9 10
> x[!duplicated(x)]
[1] 1 2 3 4 5 6 7 8 9 10
>
Question: Hello, I am new to variant calling, recently I have been ask to analyze some vcf files, most of them are using hg19.fa
but some have build37.fa in the header, it seems to be from 2014 and I am not really sure which ref genome it is, is it another name for hg19?

Answer: The correct name for the genome is GRCh37. hg19 is a nickname.