cannin/biostars_qa huggingface dataset:
Question: When using HTseq-count I get an warning at the beguining :

Warning: Read SRR1503280.2 claims to have an aligned mate which could not be found in an adjacent line.

and at the end of the sam file:

74100000 SAM alignment record pairs processed.
74200000 SAM alignment record pairs processed.
74300000 SAM alignment record pairs processed.
74400000 SAM alignment record pairs processed.

Error occured when processing SAM input (line 74468230 of file accepted_hits_filtsort_Controle.sam):
'pair_alignments' needs a sequence of paired-end alignments
[Exception type: ValueError, raised in __init__.py:603]

and my outputfiles were empty:

0 HTseq_Controle_output.txt
0 HTseq_Huntington_Output.txt

I use this code :

python ~/programs/HTSeq-0.6.1/scripts/htseq-count accepted_hits_filtsort_Controle.sam \
~/database/Mus_musculus.primary.gtf > HTseq_Controle_output.txt

I am very newbie but, I really try to look at all the other similar problems people already posted in the past and could not find a solution. Maybe it is trivial and I only need a tip. Could some one help me?

Answer: To separate out alignments from paired-end datasets and single-end datasets just use samtools:

samtools view -bf 1 foo.bam > foo.paired-end.bam
samtools view -bF 1 foo.bam > foo.single-end.bam

That should suffice. You could also just pipe that into htseq-count.
Question: Hello,

I have scRNA data. I pretend to contextualise my data with bulk RNA, i.e., merge scRNA with bulk RNA. My idea is to comprise all my cells from scRNA into a single point, and then merge that data with bulk RNA data.
I've done the first step by calculating the mean of all the cells from scRNA, resulting in a single point.
I have gathered bulk RNA, as FASTQ data, from a database. Then aligned, and extracted the counts through RSEM.

My problem is that the data don't match up. This is because scRNA uses UMI, which approaches gene counting differently from RSEM's expected or TPM counts.

I don't know how to normalise the two types of data to be able to compare them. I don't even know if it's possible.
If someone could help me, or just give me a tip, that would be great!

Thank you in advance!

Answer: Not an expert in scRNA-seq, still my two cents:

1. scRNA-seq is typically 3'-end enriched while bulk RNA-seq is fragmentated full transcript/cDNA. You would also need to ensure that the downloaded data are single-end as scRNA (at least 10x) is effectively single-end.
2. There is a strong confounding effect between your and downloaded data as different specimen, labs, kits, sequencing regimes were used. If you download a couple of independent RNA-seq datasets (bulk) for the same cell type but different kits, and perform PCA, you will see that they notably cluster by dataset rather than cell type. At least this is my experience. The effect is probably even stronger between bulk and scRNA-seq.
3. scRNA-seq is zero-inflated, while bulk RNA-seq, at normal depth, is not.
4. Bioinformatic processing is notthe same.

You better do some kind of meta-analysis like checking if prominent genes that come from your scRNA are also present in bulk, or are DEG in bulk between two conditions. I do not see a reasonable way on comparing counts between bulk and scRNA.
Question: Hello

Could you please advise any database which contains coordinates of unique loci for every isoform of alternatively spliced genes? Or any tool which can count them?

Cufflinks' output contains information about isoforms' FPKM but there are now coordinates for these loci. MATS detects different splice event but gives no coordinates.

Answer: Try [the fly genome browser @ Ensembl][1]

[1]: http://www.ensembl.org/Drosophila_melanogaster/Info/Index
Question: Hi,

Dear Community,

I have a column like this:

D309
E308
G296
T297A
P415T
P415T
V457I
V457
A214G
A214
T418
I419V
P259
P259L
L191
A190
R478
R478H
.
..
...

or in other words you can say that this column is present in a very big file as column number 19. I want only those lines in which the number matches only with the next line, that is the output should be like this:

P415T
P415T
V457I
V457
A214G
A214
T418
I419V
P259
P259L
R478
R478H

I tried this command:

cut -f19 mycolumnfile.txt | uniq -d

I got this output:

P415T

As it matches with the whole line. I want only those rows in which the number matches only.

Thanks,

Waqas.

Answer: Hi,

I would do it in python, with your list of id in the file "list" (beware, not carefully tested)

import re

with open("./list", 'r') as f1:
first = True
last_int = 0
last_line = ""
for line in f1:
if(last_int == int(re.findall("\d+", line)[0])):
if first:
print last_line
print line
first = False
else:
print line
else:
first = True
last_int = int(re.findall("\d+", line)[0])
last_line = line
Question: Hello all!

I am working with a single-cell RNA seq dataset from 10x where I have two experimental conditions. Patients pre-treatment and post-treatment condition. I want to investigate the cell-type ratios differences between treatment stages. The main problem is that the number of cells sequenced/recovered in the pre-treatment stage is so low in comparison with the post-treatment stage.

I was wondering, if exist a statistical method to "normalize" the number of cells recovered between conditions to compare the ratio in a reliable way.

Regards

Answer: You can find a detailed discussion in the "Differential Abundance" Section of the OSCA book: [https://bioconductor.org/books/release/OSCA/multi-sample-comparisons.html#differential-abundance][1], especially the section about [Composition Effects](https://bioconductor.org/books/release/OSCA/multi-sample-comparisons.html#composition-effects) might be applicable to your situation. It would be important to have replicates of the conditions (which is generally recommended anyway).

[1]: https://bioconductor.org/books/release/OSCA/multi-sample-comparisons.html#differential-abundance
Question: I am going to implement a (non)-model-organism (does that matter) genome database very soon. It will present a newly sequenced animal genome with annotation to the community. We are looking for software that will support the development of such a tool without 're-inventing the wheel'.

We will need a genome-browser for sure and a database backend for the annotation data.
We also will have use for integration of

ontologies (anatomy, developmental stage, etc.)

phenotype-genotype data, including images of phenotypes and stages

authentication/authorization

some users should be able to contribute by doing remote annotation

RNA-seq and other experimental data should be visualized as tracks

there should be private and public tracks for the genome

Blast interface to the genome

So I thought the GMOD tools will do the trick, as what we are trying to build is pretty much similar to FlyBase or WormBase, etc. Now my questions, excuse me if it is a bit vague:

How good is GMOD for this purpose and what are your experiences with it mainly as a developer? Are there good and viable alternatives to it (e.g. UCSC Genome browser + their data model)? I was planning to use G/JBROWSE and the CHADO database models, is that a good idea? Can such a system be administered by
1-1.5 persons?

We are still in the phase of summing up project requirements, but please contribute your valuable experiences, every input counts!

Answer: *cough*

Have you asked your main local end users what they prefer? I mean, of course they want everything, but they might have some preference for one or the other.

But I will also point you to an interesting survey that was done by the MaizeGDB folks when they were considering this question too: Choosing a genome browser for your organism…

And they put it all into a paper too. Choosing a genome browser for a Model Organism Database: surveying the Maize community
Question: Hello everyone! I have a general question about the importance of keeping the library fragment size consistent between samples. All my libraries have fragments ranging from 200 to 1000 bp according to the Bioanalyzer, which is fine for the Illumina paired-end sequencing. However, some of the libraries have average size of let's say 300 bp, others - 350-400 bp. I was wondering if this could be a problem by introducing some kind of bias? As far as I can see as long as the alignment rate is the same between the different samples, there shouldn't be any other hidden source of bias.
In my particular example, I am doing spiked-in RNA-seq and ChIP-seq experiments and I need to prepare libraries from sonicated genomic DNA to find out the ratio of drosophila to mouse reads in the sample to determine the consistency of the spike-in. For some reason, the sonication efficiency varies slightly from sample to sample, which results in the libraries that have slightly different fragment size distribution.
I would really appreciate if you could share with me your opinion/experience/piece of advice on this matter!

PS Crossposted on seqanswers to get more feedback

Answer: I don't have a full answer to this question but since you ask for any bit of advice, here comes my two cents :

- Variation in the sonication efficiency from sample to sample, even in the same batch, is not unusual. At least it happened to several people i know. They choose to throw away the samples that differed too much from the others before doing the library preparation, but small variations (like you have) were considered acceptable.

- A way to take fragment size into account for ChIP-seq experiments is to shift reads on each strand by half the fragment length.
Question: Hello,
I'd like to filter for sites that have a DP>5 in my normal sample and DP>10 in the tumor. My vcf has two samples columns, TUMOR NORMAL. How do I do that with a bcftools expression? I see in the documentation how to apply the filter for one or all samples in a vcf row but I couldn't find a way to specify the sample name.
Thanks a lot.

Answer: Hava a look at SnpSift filter.

With SnpSift you can build any expression that evaluates to TRUE / FALSE using the attributes of the variant and or the attributes of one or more specific samples.

http://snpeff.sourceforge.net/SnpSift.html#filter

cat my.vcf | java -jar SnpSift.jar filter "GEN[0].DP > 5 && GEN[1].DP > 10"
Question: Hi all. I need a little bit of help.
I have a list of files in a folder, they look like this:

H3K4me1/Cell-Line=S2-DRSC#Developmental-Stage=Late-Embryonic-stage#Tissue=Embryo-derived-cell-line/ChIP-chip/Rep-1//Dmel_r5.32/modENCODE_304.gff3

H3K4me2/Cell-Line=S2-DRSC#Developmental-Stage=Late-Embryonic-stage#Tissue=Embryo-derived-cell-line/ChIP-chip/Rep-1//Dmel_r5.32/modENCODE_965.gff3

H3K4me3/Cell-Line=S2-DRSC#Developmental-Stage=Late-Embryonic-stage#Tissue=Embryo-derived-cell-line/ChIP-chip/Rep-1//Dmel_r5.32/modENCODE_3761.gff3

H3K9ac/Cell-Line=S2-DRSC#Developmental-Stage=Late-Embryonic-stage#Tissue=Embryo-derived-cell-line/ChIP-chip/Rep-1//Dmel_r5.32/modENCODE_3765.gff3

I would like to change their name with sed command, so the final name would be like this:

H3K4me1_modENCODE_304.gff3

H3K4me2_modENCODE_965.gff3

H3K4me3_modENCODE_3761.gff3

H3K9ac_modENCODE_3765.gff3

I know i can play with `sed s` and the `/` but I don't manage to get it to work. I would extremely appreciate if you explain a bit the code you give as an answer.

Thank you in advance!
Best wishes,

Jordi

Answer: This is all you need in terms of an expression:

sed 's|/.*/|_|gi'

e.g.

$ echo "H3K4me1/Cell-Line=S2-DRSC#Developmental-Stage=Late-Embryonic-stage#Tissue=Embryo-derived-cell-line/ChIP-chip/Rep-1//Dmel_r5.32/modENCODE_304.gff3" | sed 's|/.*/|_|gi'
H3K4me1_modENCODE_304.gff3

Pro tip, use that regular expression with the `rename` program.

Explained:

`s` = substitute

`|` = sed delimiter so as to avoid confusion with the more traditional `/`

`/` = Match the first forward slash you find

`.*` = followed by any character, any number of times

`/` = until you meet another `/`.

`|_|` = replace all the previously matched stuff with an underscore.

`gi` = globally, and case insensitive (you don't strictly need these, I include them as force of habit).
Question: I m using deseq2 for my samples such as Wild type[untreated] , and two treated condition [with vitamin D] and Retinoic acid

and my samples are in replicates ,I want to compare the differential expression between WT vs treated[Vitamin D] and WT vs Retinoic acid .so here is my code im using but Im not sure if its doing what i want to do..

condition <- factor(c(rep("WT", 2), rep("AT", 2),rep("VD",2))))`

(coldata <- data.frame(row.names=colnames(countdata), condition))`
dds <- DESeqDataSetFromMatrix(countData=countdata, colData=coldata, design=~condition)
dds

Any suggestion or help would be highly appreciated

Answer: Looks good to me. From the dds table, you will just need to extract the contrasts now (VD vs wt and AT vs WT) using the `results()` function documented [here][1].

[1]: https://rdrr.io/bioc/DESeq2/man/results.html
Question: Hello,

I have multiple fasta sequences that are like this:

>2p__scaffold_2__5799__6580__-__778568__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>2p__scaffold_2__5799__6580__+__778569__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC
>1p__scaffold_2__11235__11438__-__830827__0.00__0.00
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTTCAATCCAGGGGATCGAGGAGATCCAAAGCAGCAGAAGCGGCTCGACGATGGTGAGG
ATTCGGGATCGGATTCAGCGCTCGTCGGGACTGG
>1p__scaffold_2__33129__34129__+__811706__0.00__0.00
GCTGGCGACGGATCTA

And I want to keep just the "> + ID" (numbers after `__+/-__` and before __0.00_0.00)

So I expect an output like this:

>778568
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGCAGCCGGA
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAA
>778569
GCTGGCGACGGATCTAGGCTCAGCGCAGAAGCAACTGAGAGTCGGCGATGAGC

I searched for it and tried this:

sed 's@.*__-__@@' input.fa > output.fa

That removed `__-__` and everything before it, including the ">" that I wanted to keep.

I also tried this to remove everything between ">" and `__-__`

sed -e 's/\>//' -e 's/\__-__.*//' input.fa > output.fa

But this removed everything after `__-__`

And this, that removed __0.00_0.00

sed 's/__0.00.*$//' input.fa > output.fa

Thank you for your help.

Answer: test=">2p__scaffold_2__5799__6580__-__778568__0.00__0.00"
echo $test | sed 's@.*__[\+-]__@>@' | sed 's@__.*@@'

Some might find:

sed 's/.*__[\+-]__/>/' | sed 's/__.*//'

to be a bit more readable.
Question: Hello,

Which tool can perform downsampling for RNAseq dataset with more than 60 million paired reads in fastq? I have used seqtk but it has issues with large memory consumptions.

Thanks

Answer: The seqtk downsampling.

# 1. Downsample a fraction of reads (useful when you know the total number):

seqtk sample read1.fa.gz 0.2 > sub1.fa
seqtk sample read2.fa.gz 0.2 > sub2.fa

It takes little memory.

# 2. Downsample a fixed number (streaming mode):

seqtk sample read1.fa 20000 > sub1.fa
seqtk sample read2.fa 20000 > sub2.fa

It uses the reservoir sampling on a stream. The peak RAM equals the size of output, independent of the input.

# 3. Downsample a fixed number (2-pass mode):

seqtk sample -2 read1.fa.gz 20000 > sub1.fa
seqtk sample -2 read2.fa.gz 20000 > sub2.fa

It reads the input twice (so twice as slow). In the first pass, it finds the sampled read indices. In the second pass, it outputs reads at the stored indices. The peak RAM is about the number of sampled reads multiplied by 24, again independent of the input. You need the latest seqtk for this 2-pass mode.
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: http://unison-db.org/

Unison is a comprehensive set of protein sequences with lots of precomputed results of many types. Unison was designed primarily for rapid mining of complex combinations of protein features. There's a simple homology-based structure viewer that can map SNPs and arbitrary regions onto structures.

The Unison source and website were released under the Academic Free License, and the data are released under the Creative Commons. There's a public PostgreSQL service (Try 'psql -h unison-db.org -U PUBLIC -d unison'). Command line tools are available. Loading and computing are automated and incremental, and coordinated by a makefile.

See http://unison-db.org/tour/ and http://unison-db.org/doc/ for examples and connection information.

The site is currently hosted at Genentech, my former employer. Unison is used internally at Genentech, but not updated as often as I used to (but there's a page with versions of source databases).

-Reece
Question: Hi all,

I'm working with a vcf (v4.1) that has incorrectly formatted deletions for some reason. The insertions are fine, but the deletions are annotated as (example):

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT
2	32474671	indel.60227	A	-	.	PASS	.	GT

Notice that the ALT is `-`, when the line should have been formatted as such (example):

#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT
2	32474670	indel.60227	GA	G	.	PASS	.	GT

I have no idea how the deletions ended up like this in the vcf, but my present plan is to parse a reference genome fasta file for these positions and manually correct all the deletion annotations, so I don't have to drop them from the vcf. What I wanted to know is if there's a tool that already does this- as it stands, I'm writing a manual parser.

Answer: I wound up just parsing the vcf in Python and calling samtools to correct the deletions in the vcf. By subtracting 1 from the alleged bp_pos, calling `samtools faidx chr:del_start-del_end` for each REF, and taking the first character of that for ALT, you can fill in the blanks. I don't believe there's a tool that corrects this problem in VCFs, because I don't think this is a common (or normal) problem. Leaving this here in case anyone ever encounters the same situation.
Question: Please let me know any tools or R packages that can create a mutation diagram showing mutations in protein domains like this figure from the MSKCC cBio Cancer Genomic Portal? Thanks in advance

Answer: Pfam provides an online tool to not only generate the domain information in JSON format, but to draw the lollipop diagram using javascript as well. They have more information [here][1].

IMHO, not as pretty as cBioPortal's but it gets you close to a solution.

**EDIT / SHAMELESS PLUG**: After seeing the data available and how easy it'd be, I made my own quick tool to fetch the data and draw the diagram for me in a style similar to cBioPortal - feel free to fork it and add features: https://github.com/pbnjay/lollipops

Example output (w/ labels per the comments)

[1]: http://pfam.xfam.org/help#tabview=tab9
Question: I have been asked to do differential and pathway analysis on some mouse RNAseq data. I was provided with FASTQs, BAMs, and gene/transcript report tables from rsem. The reads are 150 nt, paired-end from a MiSeq. rsem was run against GRCm38 with Ensembl 74 annotations. I'd like to reuse as much of the existing work as possible.

I have built some interactive volcano plots with ggvis and Shiny, and I'm starting to build a JBrowse site.

**A decent fraction** (10%?) of the genes with differential expression [Deseq2 adj p < 0.05, abs(log2fc) > 1] are **minimally annotated gene models like "GmNNN"** where the Ns are digits.

Some of those gene models have richer annotations in Ensembl 94, and some have been retired by 94. I can use BLAST to manually show that some of the reads aligning to some of these models
also align to other well-characterized genes, like protein or RNA components of ribosomes.

**Are there well-documented, web-based ways to automatically show other loci that the reads could have aligned to?**

Answer: Honestly, doing the analysis again seems to me the best. You can redo everything somewhat quickly with STAR+RSEM, or you could redo **very** quickly with Salmon or kallisto.
Question: Samtools can be used to select reads above certain mapping quality.

samtools view -h -b -q 30 aligned.bam -o above.mapQ30.bam

**But, how to select a read below certain mapping quality - all aligned reads below mapQ 30?**

I know it can be done using awk. But, the pipeline gets lengthy and time consuming when first need to convert bam to sam - separate header - use awk for mapQ below 30 - add header - sam file - convert to bam.
Really, its taking a lots of time.

Thanks,

Answer: I haven't found any straight forward method, but if someone has similar problem in the future; **below is my proposed solution:**

**Note: Check and make sure the mapped alignment is sorted.**

`samtools view -h mapped.bam -o mapped.sam`

` grep '^@' mapped.sam > RG.header` # this is important when there is header info in the sam/bam file

` grep -v '^@' mapped.sam > SAM_file.sam ` # retains the non header part of the mapped file

`awk '($5 < 40)' SAM_file.sam > Below.mapQ40.sam`

`cat RG.header Below.mapQ40.sam > Below.mapQ40_mapped.sam`

`samtools view -h -S -b Below.mapQ40_mapped.sam -o Below.mapQ40_mapped.bam`

`samtools index Below.mapQ40_mapped.bam`

**If you don't care about header info (or if its absent):**

`samtools view -h mapped.bam -o mapped.sam`

`awk '($5 < 40)' mapped.sam > Below.mapQ40_mapped.sam`

`samtools view -h -S -b Below.mapQ40_mapped.sam -o Below.mapQ40_mapped.bam`

`samtools index Below.mapQ40_mapped.bam`
Question: Hi All,

I am using 20/20 plus algorithms for the identification of driver mutations in my sample. while creating Simulated features I am facing an error.

rule simFeatures:
input: output_soma_freebayes/simulated_summary/chasm_sim_summary3.txt, output_soma_freebayes/simulated_summary/oncogene_sim3.txt, output_soma_freebayes/simulated_summary/tsg_sim3.txt
output: output_soma_freebayes/simulated_summary/simulated_features3.txt
jobid: 11
wildcards: iter=3

python `which 2020plus.py` features -s output_soma_freebayes/simulated_summary/chasm_sim_summary3.txt --tsg-test output_soma_freebayes/simulated_summary/tsg_sim3.txt -og-test output_soma_freebayes/simulated_summary/oncogene_sim3.txt -o output_soma_freebayes/simulated_summary/simulated_features3.txt
which: no 2020plus.py in (/home/ateeqanees/anaconda3/bin:/home/ateeqanees/miniconda3/envs/2020plus/bin:/home/ateeqanees/miniconda3/condabin:/home/ateeqanees/anaconda3/bin:/usr/lib64/qt-3.3/bin:/usr/local/bin:/usr/local/sbin:/usr/bin:/usr/sbin:/bin:/sbin:/home/ateeqanees/.local/bin:/home/ateeqanees/bin)
***python: can't open file 'features': [Errno 2] No such file or directory***
Error in job simFeatures while creating output file output_soma_freebayes/simulated_summary/simulated_features3.txt.
RuleException:
CalledProcessError in line 202 of /home/ateeqanees/Desktop/20_20/2020plus-1.2.2/Snakefile:
Command 'python `which 2020plus.py` features -s output_soma_freebayes/simulated_summary/chasm_sim_summary3.txt --tsg-test output_soma_freebayes/simulated_summary/tsg_sim3.txt -og-test output_soma_freebayes/simulated_summary/oncogene_sim3.txt -o output_soma_freebayes/simulated_summary/simulated_features3.txt' returned non-zero exit status 2.
File "/home/ateeqanees/Desktop/20_20/2020plus-1.2.2/Snakefile", line 202, in __rule_simFeatures
File "/home/ateeqanees/miniconda3/envs/2020plus/lib/python3.6/concurrent/futures/thread.py", line 56, in run
Will exit after finishing currently running jobs.
Exiting because a job execution failed. Look above for error message

I am not sure how to solve this I am referring to the [Tutorial Documentation from here][1].
Please let me know how to solve this issue.

[1]: https://2020plus.readthedocs.io/en/latest/tutorial.html

Answer: I think the 2020plus.py scripts is not found in your current PATH. Enter something like this into your command line or .bash_rc file:

export PATH=$PATH:/path/to/2020plus

Where /path/to/2020plus is where you downloaded 2020plus on your system. If successful you should be able to type `which 2020plus.py` into a terminal and it will print out the correct path to the 2020plus.py script. See [here][1] for more details. Hopefully this fixes your issue.

[1]: https://2020plus.readthedocs.io/en/latest/installation.html#check-your-path-variable