cannin/biostars_qa huggingface dataset:
Question: I want to save my fastq files as CRAM files to save space. However, I do want to be able to get the reads back that I have in my fastq, in case I have to reallign in a different way in the future. Below I have the code that I use to do the conversions to BAM and back again. However, when I then check the number of lines in the original fastq with the converted fastq, the original is much smaller:

original: 17460203 (lines)
converted: 147610328 (lines)

The flagstat of the BAM file gives 71038660 + 2766504 in total (QC-passed reads + QC-failed reads) reads. I don't understand why the newly made has more reads, unless the conversion from BAM to fastq also adds reads that map on multiple locations as multiple reads. Should I remove all duplicated reads and then convert?

# align with STAR
STAR --readFilesIn fastq1.fq fastq2.fq \
--genomeDir Gencode_human/ \
--outFileNamePrefix example.

# fix mates with picard
java -jar $EBROOTPICARD/FixMateInformation.jar \
INPUT=example.bam \
OUTPUT=example.fixmates.bam \
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true \
SORT_ORDER=coordinate

# convert to cram with scramble from io_lib
scramble \
-I bam \
-O cram \
-r Gencode_human/GRCh38.p5.genome.fa \
example.fixmates.bam \
example.cram

# convert back to bam
scramble \
-I cram \
-O bam \
-m \
-r Gencode_human/GRCh38.p5.genome.fa \
example.cram \
example.bam

#sort by name
samtools sort \
-n \
-o example.sorted.bam \
example.bam

#convert to fastq
samtools fastq \
-@ 4 \
-1 converted.fastq1.fq \
-2 converted.fastq2.fq \
example.sorted.bam

Answer: Hello,

I guess the problem is that you align your reads. I don't know how samtools fastq handles reads that map to multiple position. It could be that they will be also multiple times in the output.

Instead I would suggest this workflow:

**1. Create a uBam file using [FastqToSam][1]**

$ picard FastqToSam F1=fastq1.fastq.gz F2=fastq2.fastq.gz O=Sample1_unaligned.bam SM=Sample1

**2. Convert to CRAM**

$ samtools view -T genome.fa -C -o Sample1.cram Sample1_unaligned.bam

To get your fastq files back do:

$ samtools fastq -@4 -1 sample1_R1.fastq -2 sample1_R2.fastq Sample1.cram

fin swimmer

[1]: https://broadinstitute.github.io/picard/command-line-overview.html#FastqToSam
Question: Hi,

I finished mapping and variant calling on these for solid data. Now I'm looking to identify large rearrangements and transposable elements. Can anyone tell me any good tools for these? One tool I have come across now is nFuse. I wanted to know what the general pipeline for this is.

And I would like to view such arrangements as Circos plots.

Thanks for the help.

Answer: Hello Jordan.

I think SVDetect can help you.
Question: Hi,

I want to parse a SAM file based on QNAME, however the python script I designed to do so takes too long.

To visualize my issue, here are the first 5 lines of my 4.25 GB SAM file:

```
HWI-D00372:292:C6U8JANXX:2:1310:9053:40752 99 chr1_gl000191_random 623 255 51M = 737 165 CTTTGGGAGGCTGAGGCAGGTGGATCACCTGAGGTCAGGAGTTCGAGACCA /

Answer: Eek, you load the whole SAM file in to RAM as a single gigantic string, which you then regex search n times where n is the number of QNAMEs you're interested in! Thats going to need a LOT of memory and will be crayz slow! Try this:

mate_pairs = {}
QNAME_list = set(QNAME_list) # make the checking a tiny bit faster by making the list a set
with open(pathToSam, 'rb') as samFile:
for line in samFile:
QNAME = line.split('\t')[0]
if QNAME in QNAME_list:
try:
mate_pairs[QNAME].append(line)
except KeyError:
mate_pairs[QNAME] = [line]
with open(parsed_SAM_file_path, "wb") as outfile:
for read,mate in mate_pairs.values():
outfile.write(read + '\n' + mate + '\n')

This will fail if you dont have exactly two reads per QNAME (which is probably a good thing..) and should take about 5-20 minutes to run.

EDIT:

Also, dont give up on Python just because your first foray was thwarted by a pretty sweet grep!

Yes you can solve many easy problems with grep/cat/sed/awk one-liners - but I always found there are more -weird -flags in those tools to -remember, than there are programming paradigms that will never change.
Question: I've performed a nucleotide BLAST of a multi-fasta file against itself in order to find sequences which share highly similar flanking sequence, indicative of a duplication event. Given the .xml file produced from this what I need to do is determine cases where significant BLAST hits come not from the same sequence (query and subject are the exact same) but from different sequences. Essentially what I'm trying to do is parse the .xml output and get it to print details about the HSPs as well as IDs for the subject and query in each of those HSPs. I tried playing around with `f.write(header.query_id)` and `f.write(parameters.query_id)` but those just spit back an error. I suspect those need to be included in a separate for statement but right now I'm very lost.

I'm going off code based on section 7.4 of the Biopython 1.65 manual:

```
from Bio.Blast import NCBIXML
import sys

#>python screener.py inputfile.xml outputfile.txt
input_file=str(sys.argv[1])
output_file=str(sys.argv[2])
result_handle=open(input_file)
blast_records=NCBIXML.parse(result_handle)
#screens an XML file based on the specification below and prints those results to a file
with open(output_file, "w") as f:
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect< 0.01:
f.write('****Alignment****' + "\n")
f.write('sequence:'+ str(alignment.title) +"\n")
f.write('length:'+ str(alignment.length) +"\n")
f.write('e value:'+ str(hsp.expect) +"\n")
f.write(hsp.query[0:75] + '...'+"\n")
f.write(hsp.match[0:75] + '...'+"\n")
f.write(hsp.sbjct[0:75] + '...'+"\n")
```

Any suggestions would be, as always, greatly appreciated.

Answer: Ok, played around some more and solved the problem. What was needed was:

f.write(str(blast_record.query))

So the full code looks like this:

```
from Bio.Blast import NCBIXML
import sys

input_file=str(sys.argv[1])
output_file=str(sys.argv[2])
result_handle=open(input_file)
blast_records=NCBIXML.parse(result_handle)
#screens an XML file based on the specification below and prints those results to a file
with open(output_file, "w") as f:
for blast_record in blast_records:
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect< 0.01:
f.write('****Alignment****' + "\n")
f.write('sequence:'+ str(alignment.title) +"\n")
f.write('ID:' + str(blast_record.query_id) + "\n") #also added
f.write("Query Sequence: " + str(blast_record.query) + "\n")
f.write('length:'+ str(alignment.length) +"\n")
f.write('e value:'+ str(hsp.expect) +"\n")
f.write(hsp.query[0:75] + '...'+"\n")
f.write(hsp.match[0:75] + '...'+"\n")
f.write(hsp.sbjct[0:75] + '...'+"\n")
```
Question: Hello everyone,

This is the first time for me to analyze scRNA data,
compared to the example in the workflow (Seurat), I notice that I have a large amount of mtDNA,
I was thinking to incread the percentage of mtDNA so that I don't loose informations but on the other end I don't want to have bad results, do you guys have any suggestion about it ?

code:
seurat_object <- subset(seurat_object, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

my cells are leukaemia cells and here you can see the image

https://imgur.com/ZNHaMCU

Thank you

Answer: [scater](https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html#3_quality_control) has some nice QC functions that will determine which cells are outliers for QC each metric and remove them. It works on `SingleCellExperiment` objects, which are simple to convert from/to `Seurat` objects. This *feels* a bit better than just arbitrarily setting the thresholds yourself, though the results will be pretty similar.

But as previously mentioned, a threshold around 10% for mitochondrial reads is pretty standard - just make sure most of your obviously dead cells are removed.
Question: Hi,

I am using methods of analysis where it is necessary that the LD and effect sizes are computed with respect to the same reference allele, while in my case they are not since I am using an available panel to compute LD with Plink and my own data to compute effect size. To align the two measurements, I would have to compare each allele with those of the reference and flip the effect sizes in cases where the alleles are flipped, but things get messy when including both SNPs and INDELs in the comparisons. To avoid this step, I was hoping it could be possible to specify a reference allele directly when computing LD in Plink like this:

plink --bfile 1000Genomes --r2 --extract list.snps --reference-allele list.reference.allele

If this is not, are there any tools out there to easily compare alleles from two datasets and adjust the affect size of the study if the alleles are flipped compared to a reference file such as 1000 Genomes? I have previously used ['QCGWAS' - CRAN][1] which was really nice to do this, but the files I have now are all very large BED files that take too long to load into R, that is why if Plink could to this it would be much faster...

Many thanks for any suggestions and ideas!

Fra

[1]: https://cran.r-project.org/web/packages/QCGWAS/QCGWAS.pdf

Answer: From Christopher Chang

> Yes, `--a2-allele` lets you choose the reference alleles for the current run.

Based on Christopher Chang's answer, to compute LD in Plink using a specific allele as a reference we can use:

plink --bfile 1000Genomes --r2 --extract list.snps ---a2-allele list.reference.allele
Question: I have genomic positions in a tab delimited file as follows:

Chr Start End
2 10262865 11801950
2 9637403 10927601
2 25141434 27157396
2 31181368 38044662
2 127808684 129276101
2 236957807 238896574
2 237172736 238896574

**First I want to find overlapping segments within these genomic segments in the same file (edit)**. Then I wish to plot these such that the segments do not overlap on the plot. Is there a tool that can aid in::
1 Finding overlapping regions within the same file
2. Visualising the segments in a non-overlapping manner?

My attempts at plotting this in R results in the segments occupying the same position

Answer: using the tool developped for https://www.biostars.org/p/336589 http://lindenb.github.io/jvarkit/Biostar336589.html

https://gist.github.com/lindenb/08f23bb31c230878d1d744ef64f8636d
Question: I used BLAT to find the match of query B-cell receptor transcript. What does arrow+box (shown in green box) mean in UCSC genome browser? Can anyone help?

[the screenshot is here][1]

[1]: https://ibb.co/RSbdzG6

Answer: The arrow+box (which we refer to as double horizontal lines) represents an alignment with gaps in both the target genome and the query sequence. There are three ways the alignment can be drawn in BLAT results:

- **Solid black bar**: The query matches the genome entirely in that
region.
- **Single line arrows**: (as seen in the screenshot as well) These
represents alignments with gaps in the genome, but not the query. If
you were BLAT'ing mRNA, you would expect to see all the mRNA exons as
black bars and the introns as single line arrows. The arrows
represent the direction of the match.
- **Double horizontal lines**: Gaps in both the query and the target
genome. These should occur less frequently as they represent a poorer
match. There is sequence for both the query and the target genome
that differs here. The length of the double horizontal line
represents the length of the gap on the target genome, but to find
out the query length you would have to click into the alignment
itself.

If you zoom into your alignment you should see some solid alignment blocks as well. Usually these alignments are at the edges, with large double gaps in between. However, when you are zoomed out (1.3M bp), it can be hard to see these alignment blocks.

There is also some additional color coded information as part of the BLAT display. If you click on the track and scroll all the way to the bottom you should see it under the Description section.

If you have any further questions, you may reach out to the Genome Browser help desk by emailing genome@soe.ucsc.edu.
Question: I'm looking for a quick ftp-based solution to downloading FASTA canonical and splicing variants data from Uniprot:

Go to uniprot.org, type in Homo sapiens, click Download, select the format "FASTA (canonical & isoform)", click Go.

Can I automate this process for all organisms in the Uniprot database using ftp in one nice command line utility?

Answer: You can download precomputed reference proteomes by ftp. They include canonical and fasta sequences:

ftp://ftp.expasy.org/databases/uniprot/current_release/knowledgebase/reference_proteomes/

( see also http://www.uniprot.org/help/reference_proteome )

If this is not what you want, you can indeed automate the http download, as described in the UniProt help page for programmatic access, in particular Retrieving entries via queries

e.g. http://www.uniprot.org/uniprot/?query=organism:9606&format=fasta&include=yes&compress=yes

Please note that large http downloads are much more vulnerable and likely to fail than ftp downloads.

Please don't hesitate to contact the UniProt helpdesk if you have any additional questions.
Question: The package upsetr was installed with conda. It **works** in my **terminal**

conda activate upsetr
conda deactivate

I was trying to invoke conda in a **Makefile** , something like

SHELL=/bin/bash
test:
conda activate upsetr && conda deactivate

but it doesn't work:

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Currently supported shells are:

I can reproduce this error by piping the command into bash

$ echo 'conda activate upsetr && conda deactivate' | bash

CommandNotFoundError: Your shell has not been properly configured to use 'conda activate'.
To initialize your shell, run

$ conda init

Why does it fail ?

Answer: Non-interactive shells [do not run *.bashrc*](https://www.gnu.org/software/bash/manual/html_node/Bash-Startup-Files.html#Bash-Startup-Files).

So you could e.g. add the `conda init` stuff to $BASH_ENV, or alternatively make your sub-shell interactive (with all the side-effects that implies) by piping into `… | bash -i` or in the *Makefile* use `SHELL = /bin/bash -i`.

For the (GNU) *Makefile*, probably the best bet would be to avoid those side-effects via

```make
SHELL = /bin/bash
export BASH_ENV = $(HOME)/.conda.bash_env
```

Where *~/.conda.bash_env* is a one-liner: `eval "$(/path/to/bin/conda shell.bash hook)"`

Finally, you will need to ensure each command in a Make recipe that uses `conda` is actually executed via $(SHELL) rather than spawned directly by `make` via `exec(2)`/etc. The easiest way to do that is to make sure that each of those command lines contains at least one shell meta-character, e.g., just by putting a semicolon on the end.

As your `conda activate …` commands will be setting $PATH for the following commands in the recipe, this will somewhat take care of itself as you will want to ensure that both commands are run by the same shell, either using `.ONESHELL` or the usual syntax for shell script snippets:

```make
analyse:
conda activate upsetr; \
upsetr foo bar
```
Question: Hi all,

I am trying to do a couple of operations in automatic. Basically I have a directory of fastq PE files from the same sample e.g.:

L16-24MG-A_S14_L001_R1_001.fastq.gz L16-24MG-A_S14_L003_R1_001.fastq.gz
L16-24MG-A_S14_L002_R1_001.fastq.gz L16-24MG-A_S14_L004_R1_001.fastq.gz

L16-24MG-A_S14_L001_R2_001.fastq.gz L16-24MG-A_S14_L003_R2_001.fastq.gz
L16-24MG-A_S14_L002_R2_001.fastq.gz L16-24MG-A_S14_L004_R2_001.fastq.gz

I am trying to cat the fastq file together and then run bowtie so at the and having only one .bam file.
This is my script so far but it's not quite working. In fact I am not even able to obtain a combined fastq file before bowtie. Can you help me out please?

for i in $(ls *R1*.gz) do cat *R1* > ${i%.R1_combined.fastq}.gz done

for i in $(ls *R2*.gz) do cat *R2* > ${i%.R2_combined.fastq}.gz done

gunzip *.gz

for i in $(ls *.fastq | rev | cut -c 13- | rev | uniq)

do

bowtie /home/casaburi/ufrc/hybrid_pacbio_global/rsem/final_assembly_cdhit100 \
-1 ${i}_R1_combined.fastq -2 ${i}_R2_combined.fastq \
--all --best --strata -m 300 --chunkmbs 512 -S -p 10 | samtools view -F 4 -S -b -o ${i}.bam

done

Answer: Try: `for i in $(ls -1 L16*R1*); do cat ${i%%_R1*}_R1_001.fastq.gz >> ${i%%_S14*}_combined.fastq.gz; done`
Question: **Goal**

I am trying to use the FDist2 methodology. [Lositan][1] seems to not be working anymore. [Arlequin](http://cmpg.unibe.ch/software/arlequin35/Arl35Downloads.html) has an implementation of FDist2, so I am now trying to use FDist2 via Arlequin via the command line.

Note that I am on MAC OSX. If I use the command line, I can also do it on Linux.

**Issue**

I will have to specify in the .ars file that I want Fdist2 to be performed. Each specific requirement has a keyword that has to be associated to a value (typically o or 1) in the .ars file. For example `ComputeStandardDiversityIndices=1` and `LocByLocAMOVA=1` can be set to 1 as shown if we want these specific statistics to be computed. In the `arlsumstat_readme.txt`, there is a list of keywords called `Keywords in ssdefs.txt` (I don't really know what ssdefs.txt is) but these keywords don't match those used in the example files found in the .zip file for Arlecore (typically even AMOVA is absent) and none of these keywords seem to be described as what I want to do. I can't find any of these keywords in the manual either.

In the .zip file for Arlecore, there are a number of examples, incl. one called detSel which sounds exactly like what I would like to do but for some reason the .ars file is absent!

**Questions**

My questions are

- Where can I find these keywords?
- More specifically, what is the keyword for using Fdist2?
- Could you provide me with a reproducible example of a FDist2 analysis?

[1]: http://popgen.net/soft/lositan/

Answer: Facing similar issues, [Whitlock and Lotterhos 2014][1] have written their own implementation of Fdist2. It is available via [this Dryad repository][2].

It is not quite really user friendly but it comes with a R wrapper code that is nicely commented.

[1]: http://onlinelibrary.wiley.com/doi/10.1111/mec.12725/full
[2]: http://datadryad.org/resource/doi:10.5061/dryad.v8d05
Question: Hi,

Would you please recommend me a good software for searching domains in a protein or coding nucleotide sequence? (except NCBI conserved domain search). Your help is greatly appreciated.

Answer: Did you consider [PFAM][1] database?

**'Domain Organizations'** are very beautifully prepared. Look at [this][2]

You can search what domains are present in a protein sequence with `hmmscan` tool, searching the sequence against a [HMM][3] database of PFAM

[1]: http://pfam.xfam.org/
[2]: http://pfam.xfam.org/family/pf00406#tabview=tab1
[3]: ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release
Question: I have the following input:

(1) List of protein-coding gene models in GFF3 format
(2) List of protein features of these genes (e.g. protein domains or transmembrane regions), with coordinates on the protein sequence-level.

As output, I would like to have an image that shows all gene models drawn to scale side-by-side, with protein features mapped onto coding exons. Different feature types should be drawn in different colors (configurable), and features are allowed to span exon-exon junctions. Ideally, the generation of such an image can be completely automated.

Example image:

Can anyone suggest a program that can do this?

Answer: I ended up implementing two new BioPerl modules myself that, when used in combination, solve this problem. I just uploaded both modules to GitHub:

Bio::Graphics::Glyph::decorated_gene

Bio::Draw::FeatureStack

Quoting from the module descriptions:

Bio::Graphics::Glyph::decorated_gene - A GFF3-compatible gene glyph
with protein decorations.

This glyph extends the functionality of the
Bio::Graphics::Glyph::gene glyph and allows to draw
protein decorations (e.g., signal peptides, transmembrane domains,
protein domains) on top of gene models. Currently, the glyph can draw
decorations in form of colored or outlined boxes inside or around CDS
segments. Protein decorations are specified at the 'mRNA' transcript
level in protein coordinates. Protein coordinates are automatically
mapped to nucleotide coordinates by the glyph. Decorations are
allowed to span exon-exon junctions, in which case decorations are
split between exons. By default, the glyph automatically assigns
different colors to different types of protein decorations, whereas
decorations of the same type are always assigned the same color.

and

Bio::Draw::FeatureStack - BioPerl module to generate GD images of
stacked gene models

FeatureStack creates GD images of vertically stacked gene models to
facilitate visual comparison of gene structures. Compared genes can be
clusters of orthologous genes, gene family members, or any other genes
of interest. FeatureStack takes an array of BioPerl feature objects as
input, projects them onto a common coordinate space, flips features
on the negative strand (optional), left-aligns them by start
coordinates (optional), sets a fixed intron size (optional), removes
unwanted transcripts/isoforms (optional), and then draws the so transformed
features with a user-specified glyph. Output images can be generated
in SVG (scalable vectorized image) or PNG (rastered image) format.

Here is an example output of FeatureStack:
Question: I am looking for tools that would help process multiple alignment of protein sequences and extract polymorphic sites only--based on reference sequence. Well, what I have in mind is an output of something matrix-like like so:

```
Position 2 10 20 45 60 63
RefSeq Val Leu Ala Ser Phe Thr
Seq1 Met Tyr
Seq2 Gly Pro
Seq3 Ala Arg
```

Any suggestion for tools or scripts(preferably python) that can help me achieve this fast?

Thanks.

Answer: This looks like a task you can write up a script by yourself before finding a tool to do the job for you.

Assuming you have sequences saved in a python list and the reference is the last one, here's something you could do...

```
seqs = do_your_magic_to_get_the_sequences()

ref = seqs[-1]
diffs_top, diffs = [], []
all_diffs = set()

for s in seqs[:-1]:
diffs = []
for i, c in enumerate(s):
if s[i] != ref[i]:
diffs.append(i)
all_diffs.add(i)
diffs_top.append(diffs)

for d in all_diffs:
print "%-3d" % int(d+1),
print

for c in ref:
print "%-3s" % c,
print

for i, s in enumerate(seqs[:-1]):
for j, c in enumerate(s):
if j in diffs_top[i]:
print "%-3s" % c,
else:
print "%-3s" % " ",
print
```

The print out of an example data would then be like:

```
1 2 3 4 5 6 7 8 9 10 11 12 13 14
M L Y M V A R V A Y K K N P
V P R F M G L
A B K L P P K B C
N P Y F C C L C
```
Question: Hi,

I have PDX RNA-seq data.I want to separate the reads of mouse and human,before comparing with tumor RNA-seq data. Is there any tool to do that?

After removing the mouse reads I want to do differential expression between primary tumors and PDX tumors.

Thanks,
Ron

Answer: [Xenome][1] works pretty well.

[1]: http://www.ncbi.nlm.nih.gov/pubmed/?term=22689758