cannin/biostars_qa huggingface dataset:
Question: Hi all,

I have a SAM file from paired end sequencing data, I wrote a short script (below) to filter some reads (reads that map to bacteria):

awk 'FNR==NR { a[$1]; next } !($3 in a)' contigs_bacteria.txt <(samtools view -h aligned_reads.sam) > aligned_reads_no_bacteria.sam

This works fine, however when I check count of unmapped reads in the original sam file (aligned_reads.sam) vs. the filtered one, the # of unmapped reads are different:

samtools view -c -f 4 aligned_reads.sam
316971

samtools view -c -f 4 aligned_reads_no_bacteria.sam
293621

Contigs_bacteria.txt is just a text file that looks like this, and I use these values to filter the sam file, as I know reads with these values come from bacteria.

201
908
925
881
895
921
889
913
875
922
923
915
878
910

I guess it is also filtering out some unmapped reads in the process. or maybe this is due to the fact that they are paired and if only one mate meets the criteria, that pair is lost. I do not want to lose the unmapped reads as they can potentially be reassembled. Any ideas/help is appreciated.

Answer: Hello https://www.biostars.org/u/51170/ ,

if one read of a pair can be mapped and aligned but the other one not, than the one that can not be mapped is assigned to the same chromosome/sequence name like the one that can be mapped, but its flagged as unmapped. So if you exclude chromosomes/sequence names it is normal, that the number of unmapped reads decreases.

fin swimmer
Question: I am looking for a solution in both Bash/Linux and Python.

Also,

If more than 40%[Aa]’s or [Tt]’s present, remove the sequence and the description.

e.g. file.fasta

>seq1
GGCAGAGGCCCCCTAGCCCCGCCCGCGCCATGGTCAGGCACGCCCCTCCTCATCGCGGGGCACAGCCCGGCGGGTAGCCCCAGCGCTGGAGGCGGGCGGGGCCGGCCGGCGGAGGCCTGAGCAGCAGCCCAGCGCGGGCCGCCGAGACACCATGAGAGCCCCCACACTCCTCGCCCCACCGGCCCTGGCCGCACTGGGCACCGCTGGCCGGGCGGGTGGGTGCCCC

>seq2
CCACTGCACTCACCGCACCCGGCCAATTTTTTTTGTGTTTTTAGTAGAGACTAAATACCATATAGTGAACACCTAAGACGGGGGGCCTTGGATCCAGGGCGATTCAGAGGGCCCCGGTCGGAGCTGTCGGAGATTGAGCGCGCGCGGTCCCGGGATCTCCGACGAGGCCCTGGACCCCCGGGCGGCGAAGCTGCGGCGCGGCGCCCCCTGGAGGCCGCGGGACCCCTG

After processing, output should be:

>seq1
GGCAGAGGCCCCCTAGCCCCGCCCGCGCCATGGTCAGGCACGCCCCTCCTCATCGCGGGGCACAGCCCGGCGGGTAGCCCCAGCGCTGGAGGCGGGCGGGGCCGGCCGGCGGAGGCCTGAGCAGCAGCCCAGCGCGGGCCGCCGAGACACCATGAGAGCCCCCACACTCCTCGCCCCACCGGCCCTGGCCGCACTGGGCACCGCTGGCCGGGCGGGTGGGTGCCCC

seq1 is fine but seq2 and its description have to be removed from a file.

e.g.

AGAGCTAGAAGGGG - ok

AGAAAAAAAAGAGG - remove -- more than 6 A sequentially

ATTTTTTTTATGATG - remove -- more than 6 T sequentially

AGTAGTTAGGGGGG - ok

ACAGAAACAGAATG - remove -- more than 40% of A

TCTGATTTATTATTG - remove -- more than 40% of T

Are there any tools available on the market for these purposes? I would like to get a solution via tool as well use python/Linux.

Answer: **A pythonic way to achieve the same:**

from Bio import SeqIO
import getopt, sys, re

def usage():
print "Usage: remove_seqs.py -i -o "

try:
options, remainder=getopt.getopt(sys.argv[1:], 'i:o:h')

except getopt.GetoptError as err:
print str(err)
usage()
sys.exit()

for opt, arg in options:
if opt in ('-i'):
input_file=arg
if opt in ('-h'):
usage()
sys.exit()
elif opt in ('-o'):
output_file=arg

out=open(output_file, 'w')

for record in SeqIO.parse(input_file, "fasta"):
# check for 6 or more A's or T's
search_hexamers=re.search('A{6,}|T{6,}',str(record.seq))
# remove those sequences
if not search_hexamers:
# Calculate percentage of A's or T's
a = (sum(record.seq.count(x) for x in ['A'])/float(len(str(record.seq)))) * 100
t = (sum(record.seq.count(x) for x in ['T'])/float(len(str(record.seq)))) * 100
# Check of percentage of A's or T's is more than 40 % and remove
if not a > 40 or t > 40:
out.write(">"+record.id+"\n"+str(record.seq))
out.close()

**USAGE**

remove_seqs.py -i -o

**Dependency**: Biopython

**Input:**

>ok
AGAGCTAGAAGG
>remove
AGAAAAAAAAGG
>remove
ATTTTTTTTAGGA
>ok
AGTAGTTAGGGG

**Output**

>ok
AGTAGTTAGGGG

**Note**:

>ok
AGAGCTAGAAGG

above sequence was removed as the percent of A's crossed 40%
Question: Hi,

I am using DSK to count the k-mers that are present in multiple fasta files. I would like to get the list of all solid k-mers that were found in all files.

The command line interface has a "-solid-kmers-out" argument which, I assume, should output this list in a file. However, even when I set a value, it does not seem to output a file. Am I using it correctly?

Also, is there a programmatic way to go from the integer representation of k-mers to their string representation?

Thanks for your help,
Alex

Answer: Bonjour Alexandre,

To output the solid kmers in string representation, you can either use the dsk2ascii tool provided with DSK, or parse the hdf5 file using h5dump. See https://github.com/GATB/dsk#results-visualization

To integrate DSK results in a C++ program, it is currently recommended to use the [GATB-core][1] library. Have a look at [the doxygen page][2] and also the [example code][3] in the GATB-core repository that demonstrate some custom processing:

* [kmer10.cpp][4] is the most up to date version of that snippet.
* [kmer16.cpp][5] shows how to get separate kmer counts (one abundance per file) over multiple files (*multi-bank kmer counting*)
* kmer11-15 are a bit more advanced

Thanks for pointing this out by the way, the "-solid-kmers-out" argument was actually made for GATB-core graph construction (to output the graph and the solid kmers in different files), and is not used in DSK. I've added a note for the future versions of DSK.

[1]: https://github.com/GATB/gatb-core/
[2]: http://gatb-core.gforge.inria.fr/doc/api/snippets_kmer.html#snippets_kmer_snippet10
[3]: https://github.com/GATB/gatb-core/tree/master/gatb-core/examples/kmer
[4]: https://github.com/GATB/gatb-core/blob/master/gatb-core/examples/kmer/kmer10.cpp
[5]: https://github.com/GATB/gatb-core/blob/master/gatb-core/examples/kmer/kmer16.cpp
Question: Hi all,

imagine that I have a pathway P (e.g. from KEGG). The only information I have is the topology of P (nodes, edges, and direction of the edges).

Say that now I want to inhibit (modulate, if you wish, for a more general case) P. I achieve my inhibitions by formally removing edges from P. Intuitively, removing an edge that is upstream should cause a higher degree of inhibition, and possibly the same would happen for edges that depart from central nodes, and so on. However, if you take a KEGG pathway you will soon realize that the casuistics are huge, and non-trivial.

Are you guys aware of some method or approach to account for all this? As I said, I only have the pathway topology, and I want to estimate the level of pathway inhibition I would achieve by removing certain edges. I'm not an expert on the field: is this what you call "network dynamics"? Any hints will be more than welcome!

Thanks.

Answer: **Signaling Pathway Impact Analysis** (SPIA; see package here) could be an approach that could be similar to what you are looking for.

Another approach is to treat your pathway as a network and perform network level **perturbation analysis**. NEXCADE or your network library of choice with a bit of scripting and stats could do such analysis. See related manuscripts 1 and 2
Question: Hi, everyone.

I have transcript IDs as below, what does `.` mean? Are these two ID are same?

ENST00000318325.6
ENST00000318325.2

Answer: Hello,

the number after the . is the version of the transcript.

See [here][1].

When talking/writing about transcript one should always add the information about the version and the reference genome.

fin swimmer

[1]: https://www.ensembl.org/info/genome/stable_ids/index.html
Question: Hello. I have numerous bam files which I need to merge into a single bam file. How do I know whether I need to first index and sort the bam files?

Answer: You should sort your files, as `samtools merge` expects sorted bam files, either by position (default) or by name (`samtools merge -n`). You can check if you bam are already sorted by looking at the SO flag:

samtools view -H file.bam | grep "SO:"

> @HD	VN:1.0	SO:coordinate
Question: Hello,
I am having a hard time opening a very large bedgraph. I have been suggested to split my bam file by chromosome with chrom-bed.py but it didn't work.
Is there any other alternative?
Thanks,
GP.

Answer: in this other answer, Aaron Quinlan stated:

bamtools has a "split" command for exactly this purpose

I can only add that I've just tried it with this simple command

bamtools split -in file.bam -reference

and it works like a charm. the bam file gets split into different bam files, which are suffixed with .REF_xxx.bam by default, which is very convenient.
Question: According to the Variant Effect Predictor tool from Ensemble, the consequence types upstream gene and downstream gene variants are located within 5 Kb from the transcript. I would like to create 2 categories, "genic" and "non genic" and include up- and downstream gene variants (as well as all other consequence types from ENSEMBL) into one of them. Where to include up- and downstream variants?. Thanks very much

Answer: Hi there,

Ensembl uses Sequence Ontology (SO) terms tho describe the consequences of variants. The upstream/ downstream consequence type is defined as being within 5kb of the transfer:
http://www.sequenceontology.org/browser/current_svn/term/SO:0001631

So, a particular variant could be an upstream variant of all of (or a subset of) the different transcripts of a gene.

For example, if you consider the different transcripts of BRCA2 (http://www.ensembl.org/Homo_sapiens/Gene/Summary?db=core;g=ENSG00000139618;r=13:32315474-32400266) a single variant might be an upstream variant of the BRCA-001 transcript, but not of the BRCA2-003 transcript.

Technically the upstream/downstream variants fall outside of the defined gene boundaries, but the decision whether to include them in one category or another for a particular analysis really depends on the analysis you wish to perform and the definitions of variants you wish to apply during your analysis.

Best wishes

Ben
Ensembl Helpdesk
Question: Hi,

I've a bed file containing several region of interest and I want to pick randomly N positions in these regions. I know that "bedtools random" can pick random positions into the genome, but I don't find any tool to perform that on specifi region described in a bed file.

Thanks

Answer: If you want to sample `N` elements without replacement from your set of input BED elements, replace `N` with your value:

$ shuf input.bed | head -N | sort-bed - > answer.bed

Or you can use my `sample` program on Github: https://github.com/alexpreynolds/sample

$ sample -k N -p input.bed > answer.bed

This also offers sampling with replacement and some other features that might be useful, depending on the experiment.

If you want to sample `N` elements without replacement from a population of individual bases made up from your set of input BED elements, use the `bedops --chop 1` command to split your input BED file into individual bases, then do the usual `shuf` or `sample` operation:

$ bedops --chop 1 input.bed > perBaseInput.bed
$ sample -k N -p perBaseInput.bed > perBaseAnswer.bed

You could pipe the per-base output from `bedops` into `shuf` directly, but if you want to draw multiple samples, it can be a waste of time to regenerate the per-base data each time. Saving to a separate file makes for easy re-use (re-sampling).
Question: Wondering about easiest way to get a BED file of regions +/- 2kb around TSSs through table browser?

Answer: If you want to do this on the command line, you can do the following to get `hg38` refSeq annotations (adjust for your genome of interest):

$ wget -qO- http://hgdownload.cse.ucsc.edu/goldenPath/hg38/database/refGene.txt.gz | gunzip -c - | awk 'BEGIN{ OFS="\t" }{ print $3, $5, $6, $13, $10, $4 }' - | sort-bed - > refGene.bed

Then split them by strand and pad around the stranded-start position of the annotation:

$ awk '($6 == "+") { print $0 }' refGene.bed | awk 'BEGIN{ OFS="\t" }($2 > 2000){ print $1, ($2 - 2000), ($2 + 2000), $4, $5, $6 }' > refGene.tss.for.padded.bed
$ awk '($6 == "-") { print $0 }' refGene.bed | awk 'BEGIN{ OFS="\t" }($3 > 2000){ print $1, ($3 - 2000), ($3 + 2000), $4, $5, $6 }' > refGene.tss.rev.padded.bed
$ bedops --everything refGene.tss.for.padded.bed refGene.tss.rev.padded.bed > refGene.tss.padded.bed

Generally, that could give problems if a padded TSS ends up going outside chromosomal bounds. You could filter them out with a mix of Kent and BEDOPS tools:

$ fetchChromSizes hg38 | awk '{ print $1"\t0\t"$2; }' | sort-bed - > hg38.bounds.bed
$ bedops --element-of 100% refGene.tss.padded.bed hg38.bounds.bed > refGene.tss.padded.filtered.bed

For 2k of padding, this is probably not a major issue. I don't believe one would generally find many genes (and their TSSs) within 2k of the tips of chromosomes. Nonetheless, this would deal with this scenario in a general way.
Question: I have a transcriptome data of 10.8gb R1.fastq and R2.fastq each. I then compressed these pairs using `gzip R1.fastq` and `gzip R2.fast2`, and now the files are 2.2gb and 2.4gb. Is it possible for two compressed files to have different size when the uncompressed files are of same size?

Answer: A wild guess... Second-in-pair reads usually have base qualities that drops faster along the read compared to first-in-pair. This makes the quality line on each fastq record more variable (i.e. more random and less compressible) in R2 than in R1.
Question: Hi all,

I have GRCh37 fasta sequence & GTF annotation file and ERCC fasta sequence file (but I don’t have ERCC annotation file (GTF)). I want to build genome index file including ERCC sequence using STAR. please help me how to build index file with ERCC?

Thanks

Answer: Cat the fasta files together and then add an entry for each spike-in into the GTF file. Note that you could instead just map against the spike-ins and use `samtools idxstats` to get the counts for the spike-ins. The multimapping rate within the spike-ins and between the spike-ins and the human genome will be very low, so this will produce accurate enough results.
Question: I want to detect outlier present in my dataset including 205854 loci in 5 populations.
I have tried using LOSITAN on windows but unfortunately it was not working and showing "First simulation pass to determine initial neutral set" since several hours.

can anyone please suggest me how to detect outlier of this big dataset???

Thank you in advance!!

Answer: Maybe try GraphPad Prism. Alternately, you might be able to use scikit-learn in Python.
Question: I have a set of 500 100-base sequences. I want to make new sequences by randomizing this set into 500 new sequences such that the new set has the same dinucleotide distribution. So, for each sequence in my original set, I'm cutting the sequence up into mononucleotides, putting them in a hat, and drawing them out one by one, but the randomized result should have the same dinucleotide distribution. How do I do this?

The best idea I can come up with at present is to simply repeat a purely random shuffling until I find a result that is sufficiently close to the original. (closeness may be assessed with a paired t-test or, more simply, by taking the difference between the original distribution and the random one.) Is there a better way? I could slice the original sequences by twos to jump-start the process?

Background: I'm trying to reproduce a procedure described in in Weirauch et al. Nature 2012, (http://www.nature.com/nbt/journal/v31/n2/full/nbt.2486.html?WT.ec_id=NBT-201302). In trying to evaluate a TF binding prediction algorithm's success at predicting in vivo binding sites, they start with some ChIP-Seq or ChIP-exo data. From these measurements, they pick 500 loci where the TF is credibly believed to bind, and then test the algorithm's ability to distinguish the true sites from 500 control sites. This is one of the methods they use to create control sites, and they describe it thusly: "500 randomly shuffled positive sequences, where dinucleotide frequencies were maintained." This is the whole description in the methods section, leaving a few things unclear.

1. Do I shuffle the sequences individually or as a set? I mean, do I (a) slice all the sequences in the set and put them in a hat, and then test the dinuc distribution of the new set or (b) shuffle each sequence individually, as described above?
2. How precisely is the dinucleotide distribution maintained?

Answer: This is a non-trivial problem involving Eulerian walks. uShuffle claims to provide a solution:

http://digital.cs.usu.edu/~mjiang/ushuffle/
Question: I'm using PhyML to compute phylogenetic trees. In principle, there's a parallel option (via MPI), but it doesn't work for me. Instead of spending lots of time debugging MPI, I was wondering if I could run the bootstraps independently (with a user-supplied tree) on my cluster. So the real question is: Given a list of trees in Newick format, how do I calculate the bootstrap support for the original tree?

Answer: You can use an strategy similar to what Phylip does. You can generate 1000 random input files with SEQBOOT and use these files in PhyML, not in a real parallel mode, but in parallel as starting different PhyML processes in different nodes/cores at the same time (using something like this.

At the end CONSENSE will calculate a consensus tree and give you the actual bootstrap values.
Question: Hello,

I'm trying to convert a gene (DNA) into a protein (amino acids)

I tried to convert them but I have a lot more amino acids than here for example https://alphafold.ebi.ac.uk/entry/Q8W3K0

This is what I tried https://web.expasy.org/translate/

I try to find a software where I just have to enter the sequence and I will have the result

Answer: If this is a known gene then you don't need to do this manually. e.g. https://www.ncbi.nlm.nih.gov/nuccore/NM_000207.3/ is human insulin. If you look at the record you will find a link for protein sequence `NP_000198.1` in that record. Or you could simply search for protein sequence directly using the name of the gene.
Question: EDIT: Now fixed, using filters = c("chromosome_name","start","end","biotype") instead. Question still stands though if anyone knows?

Original question:
I've written some code which I think should take chromosomal coordinates (e.g. X:1-200000) and return protein coding genes (e.g. PLCXD1 and LINC00108) within that region on NCBI36.

I'm happy for it to return a gene only partially spanning the region or only genes contained entirely within that region at this stage, but I have no idea what it is actually returning.

Here's some simplified code of what I'm doing:

rm(list=ls())
library("biomaRt")
ensembl54=useMart("ENSEMBL_MART_ENSEMBL", host="may2009.archive.ensembl.org/biomart/martservice/", dataset="hsapiens_gene_ensembl")

chr.region = c("18:19092052-30289323","16:77462979-77500915","X:146419715-146584279","4:32776556-33393589","2:187947442-188506618")

entrez.ids=vector()
entrez.count=vector()
all.results=data.frame()

for (i in 1:length(chr.region)){
filterlist=list(chr.region[i],"protein_coding")

results=getBM(attributes = c("hgnc_symbol","entrezgene", "chromosome_name", "start_position", "end_position"),
filters = c("chromosomal_region","biotype"),
values = filterlist,
mart = ensembl54)

results$region = chr.region[i]
all.results=rbind(all.results,results)

ids=unique(results$entrezgene)
ids <- ids[!is.na(ids)]

entrez.ids[i]=paste(ids, sep=",", collapse=",")
entrez.count[i]=unique(length(ids))
}

write.csv(all.results, file="all_results.csv",row.names=F)

Many thanks for any advice.

Answer: Hi the issue was that chromosomal_region should be given as 18:19092052:30289323 and not 18:19092052-30289323

try this you will get the results you want

filterlist <- list("18:19092052:30289323","protein_coding")

results=getBM(attributes = c("hgnc_symbol","entrezgene", "chromosome_name", "start_position", "end_position"),filters = c("chromosomal_region","biotype"),values = filterlist, mart = ensembl54)

when you give incorrectly formatted chromosomal_region like 18:19092052-30289323 biomart takes the first value separated by ":" which in this case is 18 (chromosome) and returns all the protein coding genes on Chr.18
Question: Aloha,

I am having trouble figuring out how to remove everything after the last '_' in the sequence headers of a fasta file.

I would like this following series of headers

>ART01B_100_M7_ID100005_1
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325_189
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005_46
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

to look like this:

>ART01B_100_M7_ID100005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>PAG05A_100_M7_ID102325
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT
>KIN05B_100_M7_ALT_ID230005
TAAGAGGAGGAATTTTTCATAGAGGATTGTCTGTAGACTTAGTAATTTTTTCTCTTCATTTAGCTGGAATTTCTTCTCTT
TTAGGGGCTGTAAATTTTATTACTACAATTCTTAATTGTCGATCTTTAGGGGTTTGGTGAGATGAATTGCCCTTATTTGT

because there are various '_' in the sequence headers, this command, for example, won't work:

cat Alt_MACSE_Output.fasta | awk -F _ '/^>/ { print $1"_"$2"_"$3"_"$4 } /^[A-Z]/ {print $1}' > Alt.fasta

Can anyone help me please?

Answer: You can do this with `sed`

$ sed 's/_[^_]*$//' input.fa > output.fa

fin swimmer