cannin/biostars_qa huggingface dataset:
Question: Hi, all!

I've obtained the enriched GO BP terms using fGSEA in R as shown below.

1: GO_POSITIVE_REGULATION_OF_VIRAL_TRANSCRIPTION 0.020378457 0.08797654 0.4748196 1.5113374 13 41 TFAP4,NELFCD,POLR2H,POLR2D,NELFA,DHX9,...
2: GO_CARDIAC_CHAMBER_DEVELOPMENT 0.005847953 0.04827840 -0.3775346 -1.7146207 0 166 TEK,FOXF1,CCM2L,TGFBR3,SLIT3,TGFBR2,...
3: GO_DNA_DEPENDENT_DNA_REPLICATION_MAINTENANCE_OF_FIDELITY 0.001461988 0.04816766 0.6526836 2.0851705 0 42 EME1,DNA2,TIMELESS,BLM,TONSL,RAD51,...
4: GO_CIRCADIAN_RHYTHM 0.626463700 0.79544296 0.2392950 0.9379309 534 194 TOP2A,EZH2,TIMELESS,LGR4,NMU,SUV39H2,...
5: GO_PHOSPHATIDYLSERINE_ACYL_CHAIN_REMODELING 0.488165680 0.69185315 -0.3223920 -0.9673269 164 21 PLA2G4F,PLA2G5,LPCAT3,PLA2G2A,PLA2G1B
6: GO_SPINAL_CORD_DEVELOPMENT 0.334183673 0.55175132 0.2905675 1.0667585 261 104 MNX1,SOX4,HOXC10,SOX11,LHX5,SOX12,...

What I'm troubling is saving the result to a file.

I've tried 'write.csv' and 'write_csv', but both failed because the presence of a list (or lists) in the result that I couldn't see.

> write.csv(fgseaRes,"test.csv")
Error in write.table(fgseaRes, "test.csv", col.names = NA, sep = ",", :
unimplemented type 'list' in 'EncodeElement'
> write_csv(fgseaRes,"test.csv")
Error in stream_delim_(df, path, ..., bom = bom, quote_escape = quote_escape) :
Don't know how to handle vector of type list.

Anyone knows how to save this to txt or csv files?

Thank you!

Answer: Use fwrite from data.table package: `fwrite(fgseaRes, file="fgsea.tsv", sep="\t", sep2=c("", " ", ""))`
Question: I seem to have more reads after alignment than before.

Before alignment

` awk '{s++}END{print s/4}' reads.fastq `
> 153265

After alignment

`samtools flagstat align.bam`
> 180051 + 0 in total (QC-passed reads + QC-failed reads)

> 0 + 0 secondary

> 26786 + 0 supplementary

> 0 + 0 duplicates

> 171567 + 0 mapped (95.29% : N/A)

> [...]

I do not understand how can that be.

Can you help?

Answer: the bam contains the **supplementary** (part of the read that maps elsewhere) + **secondary** (some other probable hits for the read) alignments
Question: Hi,
I have a text file of 7500 genomic positions with their chr, start, end and want to get their nucleotide sequences. Can someone point me some tool or any thoughts on how to do it?

Thank you,

Answer: if you use python:

import pysam
genome = pysam.Fastafile(path_to_genome+'genome.fa')
sequence = genome.fetch(chr, start, end)

path_to_genome: would be any genome you have downloaded (e.g. hg19).
Question: Dear All,

I have to deal some small molecules screening for my protein-drug docking studies. Please suggest me few small molecule databases ..also the free software can use to screen them...thank you

Answer: You may want to have a look at ChEMBL (database of bioactive drug-like small molecules) and ChEBI too!
Question: Hi all,

I have two datasets A = 7 and B = 12 and the overlapping region is C =4. I want to draw a publication standard venn diagram using a R package. How can I make it?

Answer: very good and simple tutorials from the power of Google search are:

http://www.ats.ucla.edu/stat/r/faq/venn.htm

http://rstudio-pubs-static.s3.amazonaws.com/13301_6641d73cfac741a59c0a851feb99e98b.html
Question: What is the most appropriate journal for publishing this pan-genomics tool.

Also listed at omictools.com

(except Oxford's Bioinformatics)

Answer: **JANE** http://biosemantics.org/jane/suggestions.php says:

- BMC bioinformatics
- Bioinformatics
- NAR
- ....

:-)
Question: Hi folks. I need to run a de novo short-read genome assembler (on a paired-end/mate-pair library) that prefers outputting shorter but error-free contigs rather than longer contigs/scaffolds which may be mis-assembled. What assembler or what specific setting in an assembler of choice do you recommend to yield such contigs (as error-free as possible and no contig overlappings)?

Answer: According to [this paper][1] in BMC Bioinformatics journal:

- For short read libraries (e.g. Illumina MiSeq): **CLC bio assembler** ([CLC Assembly Cell][2]) (commerical, free 2-week trial)
- For Roche 454 read libraries: **Newbler** (Roche)

These assemblers tend to break reads and contigs at repeat boundaries and place repeated elements into separate contigs. Hence we might have more conservative and better quality (less likely to be mis-assembled) contigs.

[1]: http://www.biomedcentral.com/1471-2105/15/211/abstract
[2]: http://www.clcbio.com/products/clc-assembly-cell/
Question: Hello,

I'm trying to call SNPs from a .SAM file.

I've checked the output format of .vcf files but I couldn't find what I need.

I need to extract the counts of alleles for each SNP ID.

Something like this :

rs99999999 A:123, T:0, C:345, G:0

I'm a beginner in the field so sorry if my question looks simple but I cannot find the answer.

Can anyone tell me how to do it ?

Thanks in advance.

Answer: You won't find this kind of count directly into a VCF file. But there is the allelic depth for the reference and the alternative alleles.

For example for this SNP:

chr_1 21682 . T C 150.0 . AC=1;AF=1.00;AN=1;DP=4;FS=0.000;MLEAC=1;MLEAF=1.00;MQ=56.44;QD=31.06;SOR=3.258 GT:**AD**:DP:GQ:PL 1:**0,4**:4:99:180,0

The reference allele is T and the alternate allele is C. And if you look at the AD (allelic depth) field, you will find that there is 0 reads supporting the reference, and 4 reads supporting the alternate allele.

You can use a tool like SnpSift extractFields to get this field.
Question: I wonder why I can't find any database for enzymes composed by multiple proteins as separate entities.

I'll explain better: the OST complex is a multiple-units enzyme composed by 6-7 proteins, like OST4, STT3A or STT3B, RPN1, DAD1, etc...
If I look for any of these proteins in Uniprot, I get a nice entry for each: for example, if you look at STT3A you get a lot of informations, including a description that explains that this protein partecipates to the OST complex.

However, if you look for 'OST complex' in uniprot, there is not a separate entity for this multiple-units enzyme... I wonder if you know of a database where complex enzymes are described, or if you have faced this problem before.

Answer: Does this one help? CORUM – the Comprehensive Resource of Mammalian protein complexes.

Just tried a quick search using "DAD1" and it seemed to return useful results.
Question: Dear all,

is it possible to visualize insertions in a sequence?

I have prepared a simulated sequence of the mitochondrial genome from the release hg38 by placing non human sequences right in the middle of it (position 8284). I then aligned the simulated genome to the mitochondrial index and the visualized the alignment with the integrated genome viewer IGV. However, I don't see any sign of insertions in the figure.

![enter image description here][1]

Is there a way to highlight the insertion point? Maybe by showing only clipped reads or the reads that map only on one mate?

Thank you.

[1]: http://u.cubeupload.com/gigux/Msplitpoint.png

Answer: I see evidence of the "transgene" insertion: all those **identical soft-clipped bases** centered at the position you inserted the non-human sequence. Pay attention: 1) all reads are soft-clipped at the same reference position, 2) as far as I can tell, all soft-clipped bases are identical between different reads.

Look at the picture below. The big red arrow indicates the insertion point, and the darkened rectangles indicate the inserted sequence (which I was able to determine as parvovirus by blasting them, even before you told us it was parvovirus).

However, keep in mind this visual inspection works well because you have a simple, small and with no duplications reference genome, and a simple and small insertion, without other copies of it throughout the reference genome. As WouterDeCoster [pointed above][1], there are better methods to identify structural variation events in more complex scenarios.

[1]: https://www.biostars.org/p/331663/#331694
Question: I have a large metagenomic RNA-seq dataset that I am trying to assemble to find viral sequences but it is too large for my hardware (52gb RAM). I can see that there is a lot of bacterial contamination from many different species when I BLAST reads. I want to filter out all bacterial reads so that I can assemble. Ideas?

1. Download all bacterial genomes from Refseq and try to bowtie to that (will take a long time). As well, since when has the compressed Refseq bacterial fna files reached 72gb (when combined)?!? The last all.bacteria.gz file in Refseq archive from 2015 is 2.7gb...

2. Somehow condense all bacterial genomes into non-redundant, then align?

3. Other ideas?

Answer: You cannot determine whether there is bacterial contamination simply from BLASTing reads. Similarly, it is impossible to trivially filter out bacterial contamination by mapping to all known bacteria, because viruses tend to share sequence with their hosts, and there's no guarantee that your bacteria are in the reference dataset.

You need a completely different approach. Perhaps you should assemble the data, annotate the assemblies, and then pull out contigs with genes known to occur only in viruses...
Question: Hi,

I'm having trouble removing duplicates using Picard tools on SOLiD data. I get a regex not matching error.

The reads have the following names:

22_758_632_F3

604_1497_576

124_1189_1519_F5

358_1875_702_F5-DNA

And I don't think Picard tools is able to pick these read names with its default regex.

I tried to change the default regex. This time it does not throw an error, but it takes too long and times out (out of memory). I suspect I'm not giving the right regex. Here is my command:

java -jar $PICARD_TOOLS_HOME/MarkDuplicates.jar I=$FILE O=$BAMs/MarkDuplicates/$SAMPLE.MD.bam M=$BAMs/MarkDuplicates/$SAMPLE.metrics READ_NAME_REGEX="([0-9]+)_([0-9]+)_([0-9]+).*"

Any help is appreciated. Thanks!

Answer: I was able to fix the issue, by adding -Xmx16g and increasing the RAM size. Apparently the RAM was not sufficient.
Question: Hi every one i am doing some plamid genome assembly with spades. After assembly i used SSPACE for scaffolding. But there are some gap each of the draft genome. I can fill the gap by PCR. But i want to reduce the number of gap insilico?? Can anyone suggest how to reduce the gap ?? If i do mapping the reads with contigs will it be give any promising result???

Answer: GapCloser from SOAPdenono works perfect for me.

To use it, you will need paired-end or mate-pair reads.
Question: Hello,

I wish to filter from a big FASTA file only sequences related to chrs 1-22, X and Y.
An example of FASTA sequence is:

>ENSG00000119314|ENST00000210227|PTBP3|9|-1|115024785
GGGTGGCAGGTGCCTGTAATCCCAGCTACTCCAGAGGCTGAGGCAGGGGAATTGCTTGAG
CCTGGGAGGCAGAGGTTGCAGTGAGCCGAGATTGTGCCACTGCACTCCAGCCTGGAGTCT
CACTTTGTCACACAGGGTGGAGTGCAGTGGTGTGATCTCGGCTCACTGCAACCTCTGCTT
ACCGGGTTGAGATTCTCCTGTCTCAACCTCCTGAGTAGCTGGGATTACAGGCGTGCACCA
CCAAGCCAGACTAATTTTCCTATTTTTAGTAGAGATGGGGGGTTTCACCATGTTGGCCAG...

>ENSG00000236011|ENST00000211377|GPANK1|HSCHR6_MHC_COX|-1|31616421
CCCTATTCCTACCTAACCTCCCCTCAGGACTCAGGCTCCAATGTGTTGAGCCCCAACTCC
TTCCCATAAGACTGCCACACGGTGCTTTCCTTTCCCTTCTTCAACACTCACCAATGGGAA
GCATTGGCTGGTTCTCACAGTACACACGAGGACAGTAACCAAAGTCTCCTTGCTGGTACT
TTTCCAACTGAGGTGAATACAATGGAAGGGGTTGGCAGGTAGATGTAAAGAAGAGGCAAC
TCCCTTCGCAGCCCAACCCATACCACTCTGTCCCCCACTCCTCCCACCTCTGTCCAGAGG
CCCCTTCTCTGGACTAGACGGGCTCTCAAACTTCTGTGTTGCCTTTCTTCCAATTAGGCA
GGCTACAAACCATCAGAGCCATTTGTTGTTTGTTCCTTGAGGAAGAGGCAGTCTATCACA
ACTCTCTGATTCAAGGTCTGTCTCCCTCCCTGAAAACAATCCCTTCAGGATGACCCCCAA...

>ENSG00000087494|ENST00000201015|PTHLH|12|-1|28115255
TCCGCTCACGGGCCCCGAGACCCCCGAAGTTCCCATGGAGCCTAAGATCCCCAGGAGCCA
AGCCTGCCCCGTCCCTGCGGATCAGCTTCCTAATGGGCGACCCAAGTCTATCGCAGGCGG
TGGGGATGAGGACGCTGGGTGGGAGGAGGGGAGGGGAGGCTGAAAAAGATCATCCCCCTT
GCCCTAAGGCCTCTCCCAAGACCCTGGACCCCTGCCCTAAGAGACTCAGGCCTCCCTTGC
TGCAGTGGGAGCGCAAACACCAGGGCAGGAGACTCCAGAGAAGGAGCGCATAACTCAACG
TTTGCTCTCCTGAAGCCTTATTTCTGATAAAAATTACAGAAAAGTTAGGCAGGATCCAAA
GACACCGTAATGACCAGCTCAAAGCCAAACAGACAGGACATCCAGTGCGGGTGTCTGGAT...

As you can see in the forth place after "|" there is the chromosome name:
in the first one is: 9
the second one is: HSCHR6MHCCOX
and the third one is: 12

I want to create a new FASTA file containing only 9, 12 sequences and in more generally 1-22, X and Y sequences.

What is fastest way to do it?

Thanks,

Tom.

Answer: You can use samtools for that:

#index a genome
samtools faidx human_genome.fa
#select chromosomes or regions
samtools faidx human_genome.fa chr1 chr2:1:2000 [...] > human_selected.fa
Question: Hi all,

I've come across a problem in PLINK when trying to do a Fishers exact test. The command I'm using is as follows:

plink --file test --fisher --allow-no-sex --1

And the error I get is:

```
ERROR: Locus 1:54208 has >2 alleles
Individual Ind3 Ind3 has genotype [ G G ] but we've already seen [ A ] and [ T ]
```

I've checked my file rigorously and the data is indeed 'GG' with no A's or T's nearby! I also have no missing data. The length of each line (i.e. for each individual) is consistent throughout. I've tried both tab- and space-demilited files, but no difference. I haven't found any special characters etc. either (using vi :set list).

Interestingly, I've taken Ind3 out of the file and re-run the test, but the same error is thrown up (but now obviously on Ind4, which is now on line 3).

Any ideas?

Answer: Plink requires that sites be biallelic. If ANY other individual has a nucleotide/nucleotides that make it multiallelic at that site, then plink fails.

Barring this, your file is formatted incorrectly. From the plink manual:

> Genotypes (column 7 onwards) should also be white-space delimited; they can be any character (e.g. 1,2,3,4 or A,C,G,T or anything else) except 0 which is, by default, the missing genotype character. All markers should be biallelic. All SNPs (whether haploid or not) must have two alleles specified. Either Both alleles should be missing (i.e. 0) or neither. No header row should be given. For example, here are two individuals typed for 3 SNPs (one row = one person):
>
> FAM001 1 0 0 1 2 A A G G A C
> FAM001 2 0 0 1 2 A A A G 0 0
> ...
>
> The default missing genotype character can be changed with the --missing-genotype option, for example:
>
> plink --file mydata --missing-genotype N
Question: Dear all,

"TO TRIM OR NOT TO TRIM?"

My PE RNAseq library prep of human brain tissue was made with TruSeq Illumina kit A using index 5, and I've got a few yellow warnings that I'd like to know what you'd do.

I found a yellow warning for overrepresented sequences - none are Illumina adapters/index. When I align the multiple overrepresented sequences, they mostly overlap, and when I blast that sequence this is the result:

[Homo sapiens uncharacterized LOC105378179 (LOC105378179), transcript variant X2, ncRNA][1]

Should I REMOVE this sequence using Trimmomatic, since it's overly reprersented?

For last, I have a warning on per sequence GC content.

I thought of doing a "mild" trimming of the reads using trimmomatic (`LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36`) to remove low-quality bases, and that's all. What would you recommend?

[1]: http://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_767982462

Answer: For what it's worth, I always clean my raw data, at the very least for poor-quality base calls or poor-quality reads in general. So do my colleagues. The real (and somewhat longer) answer to your question, however, is rooted in what type of data you're generating and what you want to ultimately do with it.

If you have shotgun libraries and are looking to assemble a whole genome, keeping in low-quality reads and bases increases complexity and can dramatically increase run time. It's pretty striking; I've seen assembler get 'hung up' trying to sort out the k-mer graph, and the problem can disappear once poor reads are removed. This also applies to duplicates.

If you're just mapping to a reference and calling variants, it's less of a deal nowadays than it was a few years ago. BWA's MEM algorithm can soft-clip reads to improve mapping quality, and this is useful if you have residual adapters or low-quality spans at the beginning. I see this as a secondary bonus of sorts, but I would still trim my reads.

Also, you might have a non-random distribution of k-mers or subsequences represented just based on your library prep. Imagine that you PCR a single locus and then make a library out of it. You will definitely have an overrepresentation. Scaling up, say you targeted and sequenced the exome. Again, your distribution of k-mers might be non-random because you might expect to see certain motifs overrepresented (start/stop codons, for example).

So, I would clean my raw data using a series of best practices (remove low-quality bases/reads/adapters, identify overlaps, dedup - but the dedup doesn't apply to your expression data). I would also be a bit leery to just chop bases off for no reason other than a summary report suggests overrepresentation. The question to ask is, "Will this affect the biological interpretation of my data systematically?"

I would love to hear others' thoughts.
Question: I'm just trying to use samtools mpileup and it is aborting with out showing any error for some bacterial genomes, with the file given there are two variants and the program dies.

The error is just:

```
[mpileup] 1 samples in 1 input files
Set max per-file depth to 8000
Aborted (core dumped)
```

The reference file (indexed with bwa) is: https://drive.google.com/file/d/0BzPjZ1hM-XWPMmxucnpDUjdpUFE/view?usp=sharing
The sorted bam file is: https://drive.google.com/file/d/0BzPjZ1hM-XWPUWw4NkdhM0E3Vm8/view?usp=sharing

The call is:

samtools mpileup -uf SaureusTCH1516_nt_genome.fasta 446.srt.bam

To obtained the sorted bam file:

First I trimmed pair end reads by quality and if adapters were found.

Second the trimmed reads are splitted in two different files, one for each end of the sequencing data, resulting in files without quality information.

Third I aligned the reads files to the reference with bwa aln and bwa sampe

Fourth I created the sam file with samtools view | sort

And finally I indexed the sorted bam file resulting in the file given in the link.

The version of samtools I'm using is 1.1 using htslib 1.1

Any hints on what I could be doing wrong, or what may be happening?

Thanks in advance

Answer: You have a corrupted alignment (for read #1 of HISEQ:185:HBB8WADXX:2:1105:8685:78337). I haven't determined what's actually wrong with its encoding (well, it's with the auxiliary tags), but you can simply avoid the problem with

samtools mpileup -Ruf SaureusTCH1516_nt_genome.fasta 446.srt.bam
Question: I am unable to view my data which I was able to view and manipulate with similar commands prior to this. Running on a Linux server. This is what I run and the errors I get:

```
​tabix sequence.snps.vcf.gz chr37>chr37.vcf
bgzip -c chr37.vcf > chr37.vcf.gz
tabix -p vcf chr37.vcf.gz >chr37.vcf.gz.tbi
```

ERROR:

tabix the index file already exists. please use '-f- to overwrite

I saw that I had the tbi file so thought I could just run this to see some of the file

```
tabix chr37.vcf.gz chr37:1-20000
[ti_index_load] wrong magic number
[ti_index_load] fail to load the index: chr.37.vcf.gz.tbi
[tabix] failed to load the index file
```

I removed the index file and remade it didn't fix the problem. I tried

```
tabix -p -f vcf chr37.vcf.gz > chr37.vcf.gz.tbi
[main] unrecognized preset
```

I tried moving the `-f` in front of the `-p`, removing the `-p`, and the `vcf` still get same error. I even reloaded the `sequence.snp` file and remade the files. I have no idea what the issue is. Any help would be greatly appreciated.

Answer: I believe the command you want for that last part is just:

tabix -p -f vcf chr37.vcf.gz

It automatically creates the appropriate file. When you add the redirect after it, it immediately creates an empty file and prepares to write stdout to that file. Then, when tabix checks, it sees that a file by that name exists.
Question: I have some RNA sequencing reads to align to the human reference genome. I found the genome FASTA files on both [GENCODE][1] and [ENSEMBL][2]: `GRCh38.p13.genome.fa.gz` and `Homo_sapiens.GRCh38.dna.toplevel.fa.gz`
But after unzipping them, I found that they are 3.1G and 60G respectively. Why is that? And which one should I use? (considering the purpose of the project is to detect gene fusion from the sequencing reads).

[1]: http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/
[2]: http://ftp.ensembl.org/pub/release-104/fasta/homo_sapiens/dna/

Answer: `toplevel` file from Ensembl includes haplotypes with full length of chromosome padded out using N's. That is the reason it is huge compared to GENCODE file. Use the Ensembl primary file which is equivalent to GENCODE.

From README at Ensembl:

---------
TOPLEVEL
---------
These files contains all sequence regions flagged as toplevel in an Ensembl
schema. This includes chromsomes, regions not assembled into chromosomes and
N padded haplotype/patch regions.
Question: Hi,

I am working in R environment and I have a long list of gene names which I am showing few first objects of it here:

```
[1] SPNCRNA.1436,omh5,snR95
[2] snR46
[3] snR10
[4] SPNCRNA.1651,SPNCRNA.515
[5] snR42
[6] SPNCRNA.1094,SPNCRNA.1095,SPRRNA.47,SPRRNA.48
[7] snR88
[8] SPNCRNA.497
[9] SPSNORNA.54
[10] snoR39b
```

I am wondering if there is any way to split the indexes with several gene names in them into individual ones? The list is so long.

Thanks!

Answer: to elaborate a little on RamRS answer:

```
x <- c("SPNCRNA.1436,omh5,snR95","snR46", "snR10", "SPNCRNA.1651,SPNCRNA.515")
results <- c()
for (i in 1:length(x)){
n <- 1
xi <- strsplit(x[i], ",")
results <- c(results, xi[[1]][n])
# print(xi[[1]][n], sep="")
while (!is.na(xi[[1]][n+1])){
n <- n+1
# print(xi[[1]][n], sep="")
results <- c(results, xi[[1]][n])
}
}

results
[1] "SPNCRNA.1436" "omh5" "snR95" "snR46" "snR10" "SPNCRNA.1651" "SPNCRNA.515"
```

--added some edits--

RamRS already gave you a perfect code, I still modified mine above

with this you can then save the results in whichever way you prefer.

`x` is your initial vector with the names, if it's part of a column then `x` will be `myMatrix[,(# col with x)]` or `myDataFrame$x`