cannin/biostars_qa huggingface dataset:
Question: Hi All,

I'm new to RNAseq, so I'm confused with a couple of terms.

As I understand, raw counts are generated by sequence aligning. Is it right?

And then, expected counts are generated by running RSEM with raw counts, right?

This is to correct sequence reads mapped to more than one spot or alternatively spliced variants. Is this right?

Then, expected counts are not normalized by length of genes. Did I understand correct?

I guess, this is pretty simple things, but I'm still a little confused when someone asks me..-.-v

Please correct, if I have misunderstanding.

Thanks in advance.
MH

Answer: Raw counts are what you get when you align to the genome and count how many reads (generally discounting multimappers) align to each gene.

What you refer to as "expected counts" are better called "estimated counts". This involves using expectation-maximization (or similar) to try to apportion multimapping reads in the most likely manner. This can also be used to get isoform-level metrics. Aside from RSEM, salmon and kallisto can also do this.

Estimated counts are typically presented in TPMs and/or estimated counts. TPMs are length normalized, while estimated counts are not.
Question: Hello,

I'm new here and and I need help please!

I need to download a list of all human genes with their respective gene symbol | chromosome | strand | transcription start site | Txen | and **Ensembl gene name**,

Actually I'm using UCSC table and I get something like this:

```
#hg38.knownGene.name hg38.knownGene.chrom hg38.knownGene.strand hg38.knownGene.txStart hg38.knownGene.txEnd hg38.kgXref.geneSymbol
uc001aaa.3 chr1 + 11873 14409 DDX11L1
```

And I want to know if it's possible to include ensembl gene symbol with UCSC table or with another method

Thank you in advance

Cherif

Answer: Here's an easy way to do it from UCSC's table browser:

1. In the table browser, select Ensembl Genes as your track
2. Under output format, choose "selected fields from primary and related tables"
3. Add your output file name, if you want one (otherwise, it will print to the browser)
4. Click "get output"
5. On the next page, you will get to choose your fields.
6. Under linked tables, check `ensemblToGeneName`, then press "Allow selection from checked tables"
7. The page will refresh, and you should have a new table called `hg19.ensemblToGeneName`
8. Check of name, chrom, strand, txStart, txEnd, and name2 in hg19.ensGene (or any fields you'd like)
9. In `hg19.ensemblToGeneName`, check "value", which has the description "alternate gene name"
10. Press "get output"

If you did it right, you should get a table that looks a bit like this (I took this chunk from chr1:100,000,000-150,000,000):

```
#hg19.ensGene.name	hg19.ensGene.chrom	hg19.ensGene.strand	hg19.ensGene.txStart	hg19.ensGene.txEnd	hg19.ensGene.name2	hg19.ensemblToGeneName.value
ENST00000263174	chr1	+	100111498	100160097	ENSG00000099260	PALMD
ENST00000605497	chr1	+	100111748	100155633	ENSG00000099260	PALMD
ENST00000605613	chr1	+	100133135	100135379	ENSG00000099260	PALMD
ENST00000496843	chr1	+	100148821	100160097	ENSG00000099260	PALMD
ENST00000434734	chr1	+	100163797	100164734	ENSG00000223656	HMGB3P10
```
Question: I am on the hunt for a maximum parsimony tree-building program with bootstrapping for Linux that one can parse options to easily from the command-line. Does anyone have any recommendations?

In brief, the tree-builder is to be run as a part of a shell script based pipeline that produces consensus trees for a given alignment under various conditions and finds the configuration that produces the consensus tree with the highest mean branch support value (thus the bootstrapping). My alignment has ~100 sequences of ~16kb, though these are quite highly conserved (parsimony informative sites are only a few hundred). As the earlier steps of the pipeline produce alignments in .fasta format it would be preferable that the tree-builder can read these.

I have been looking at TNT so far but find its interface rather unwieldy. The fact that it opens an interpreter when the program is called is sub-optimal as I would prefer that I can simply run the program with options from the command-line from the shell script (i.e. user$./program [options] [file]). This being said, if anyone has figured out how to use TNT like this then I would be grateful for your advice.

The other option I have looked at is PAUP* but I would prefer not to purchase this if there is a simpler program freely available.

Thank you for any help you can offer Biostars folks.

Answer: Zev referred me to this question, which provided me the impetus for registering on Biostars.

As he mentioned, PAUP* is now open source. However, I cannot find where Dave Swofford was posting the change log. The worst part is that I know I've been there before - I recall him saying that Sinauer will have the license for the GUI but the command-line version will be free, which is what you want. I posted a comment to Zev's response because I was still looking for the URL; I'll post it if I can find it (perhaps it's down right now?).

Anyway, you may consider R. Klaus Schliep's phangorn library (http://cran.r-project.org/web/packages/phangorn/phangorn.pdf) will estimate phylogenies using MP as an optimality criterion. If you need a complete command-line interface for some other stuff in a pipeline, I recommend writing a script and calling it using

Rscript --vanilla your_script.R [arguments]

When called this way, R, like other languages, will accept command-line arguments.

Hope this helps!

[Edit: typos]
Question: I am using esearch query as `$query = "SS1G_03709+AND+gene[filter]";`, but it gives me all (gene+ mRNA+genome sequences). What filter do I need to use so I only get gene sequences in my search? I tried a few filters from [here][1], but couldn't find anything to limit my search for genes.

[1]: https://www.ncbi.nlm.nih.gov/books/NBK49540/

Answer: If you are interested in just the FASTA sequences of mRNAs, you can use `elink` as follows:

esearch -db gene -query 'SS1G_01676' | elink -db gene -target nuccore -name gene_nuccore_refseqrna | efetch -format fasta

If you happen to know the unique gene ID (for example, it is 5493342 for SS1G_01676), you can skip the esearch command entirely and go directly to the elink step as follows:

elink -db gene -target nuccore -name gene_nuccore_refseqrna -id 5493342 | efetch -format fasta

Another tool that you may find useful is `efilter` which you can use to filter results based on source (refseq, genbank, etc) and/or molecule type (genomic, mrna, etc) as follows:

esearch -db nuccore -query "SS1G_01676 [GENE]" | efilter -molecule mrna | efetch -format fasta
Question: Dear all,

I have a problem with Biopython.

I give it a fasta-sequence and need to make either reversed fasta sequence out of it in a separate output file, or just complement fasta-sequence in a separate output file.

The four lines below were taken from Biopython cookbook, and the script works perfectly well.

from Bio import SeqIO

records = (rec.reverse_complement(id="rc_"+rec.id, description = "reverse complement") \

for rec in SeqIO.parse("example.fasta", "fasta") if len(rec)<700)

SeqIO.write(records, "rev_comp.fasta", "fasta")

18

But my goal is fasta-file with just complement to the initial sequence, or reverse sequence to the initial fasta sequence.

So far I have not found any way to get them in Biopython.

I can open final rev_comp.fasta, copy the sequence itself, do something like that:

my_seq[::-1] and print the result to a new fasta-file, it will give me just complement sequence.

But it won't be a pure Biopython application, and a little bit more complicated with reverse sequence only.

Is there any simple way to do it in Biopython only? Actually I don't need reverse complement sequence at all,
just either reverse or complement fasta sequences for my initial fasta-sequence.

Thank you very much!

Natasha

Answer: Here you go:

**USAGE: python complement.py -i input_fasta -o output_fasta**

from Bio import SeqIO
import getopt, sys

def usage():
print "Usage: complement.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"):
out.write( ">"+record.id+"\n"+str(record.seq.complement())+"\n" )
out.close()

Let me know if it works & don't forget to upvote ;)
Question: Hello,

I downloaded the files from: https://bigd.big.ac.cn/gsa/

The file is ended with: " *._f1.fq.gz" and " *._r2.fa.gz".

1. Is it single-end or paired-end sequencing?

2. If it is paired-end sequencing, the file should be: " *._r1.fq.gz" and " *._r2.fa.gz", not " *._f1.fq.gz".

Thanks in advance for any help!

Best,

Yue

Answer: > If it is paired-end sequencing, the file should be: "._r1.fq.gz" and "._r2.fa.gz", not "*._f1.fq.gz".

There are no hard rules regarding the labeling of paired files, which is most likely what yours are. And if there are rules, they are not followed by everyone. In your case, these files are likely `forward` (r1) and `reverse` (r2). It is not difficult to verify this after you unpack the files and type:

head *_??.fq

If both files have similar headers except where one of them has `1` the other has `2`, they are paired-end files.
Question: Good afternoon,

I'm trying to filter a VCF file that has the following dummy flag values:

- `PASS`: All filters passed
- `Fa`: Failed filter a
- `Fb`: Failed filter b
- `Fc`: Failed filter c
- `Fd`: Faield filter d

Variants can fail one or more filters. Variants that fail multiple filters will be annotated with the corresponding flags separated by semi-colon. Thus, the filter column can have one of the 5 above values, or any number of `F*` values separated by `;`.

I'd like to select all variants that either `PASS`ed or only failed filter a. How can I do this in bcftools? The `-f` option skips location that does not contain one of the listed filters, so it keeps locations that contain any of the listed filters. When I use

bcftools view -f PASS,Fa ...

I get rows that failed filter a along with other filters also. That is, the above expression matches both `Fa` and `Fa;Fb`. I tried excluding the delimiter, but that didn't work:

bcftools view -f 'PASS,Fa,;' ... #didn't work

Does anyone know how to exclude or include *exactly* a list of filters? Nothing in the `-i` or `-e` EXPRESSIONS is useful either.

This is what I'm using right now, which is `awk` mocking `bcftools`:

zcat vcf_file.vcf.gz | awk -F"\t" -vOFS="\t" '$0 ~ /^#/ {print} $7=="PASS" || $7=="Fa" {print}'

Answer: This feature was lacking in bcftools, and the developer has now fixed that with this commit: https://github.com/samtools/bcftools/commit/fea8773196878481399183fd9f711685d41e6cf9

Starting bcftools v1.10.3, it should be possible to do this sort of exact filtering.
Question: Hello All,

I have one eight lane HiSeq-4000 flow cell data.I want to demultiplex only lane 5,6,7,8 and skip other lanes.I made a sample sheet with only lane 5,6,7 & 8 and processed the bcl2fastq by the following command. But it generate lane 5,6,7 & 8 lane data along with Undetermined fastq file of lane 1,2,3 & 4.It is time more time & space consuming process.Please let me know that Is there any way to skip the unwanted line? So we can save time and space while doing bcl2fast.

Command Used:

bcl2fastq --input-dir 151112_K00155_0019_BH573GBBXX/Data/Intensities/BaseCalls/ --runfolder-dir 151112_K00155_0019_BH573GBBXX/ --sample-sheet Samplesheet.csv --output-dir unaligned --use-bases-mask Y*,I8n,Y* --barcode-mismatches 1 -l DEBUG &

Answer: Are you using bcl2fastq2? If so, I think you could use the `--tiles` parameter:

--tiles s_[5-8]

**edit**: I've updated my answer according to Devon Ryan comment bellow, previously it was `--tiles s[5-8]`
Question: hi,

please help me with this obstacle :(

with this syntax I tried to map my reads on my indexed genome

[izadi@lbox200 bin]$ **mapper.pl clean.fq -e -p Aspergillus_fumigatus -s reads_collapsed.fa -t reads_collapsed_vs_genome.arf -v**

but

parsing fastq to fasta format

discarding short reads

mapping reads to genome index

trimming unmapped nts in the 3' ends

Log file for this run is in mapper_logs and called mapper.log_2598

Mapping statistics

#desc	total	mapped	unmapped	%mapped	%unmapped

Use of uninitialized value $count2 in subtraction (-) at ./mapper.pl line 708, line 2451201.

Use of uninitialized value $count in subtraction (-) at ./mapper.pl line 708, line 2451201.

total: Use of uninitialized value $count in print at ./mapper.pl line 708, line 2451201.

Use of uninitialized value $count2 in print at ./mapper.pl line 708, line 2451201.

0	Use of uninitialized value $count in division (/) at ./mapper.pl line 709, line 2451201.

Use of uninitialized value $count2 in division (/) at ./mapper.pl line 709, line 2451201.

Illegal division by zero at ./mapper.pl line 709, line 2451201.

[izadi@lbox200 bin]$

Answer: [izadi@lbox200 bin]$ mapper.pl clean.fq -e -h -m -v -p Aspergillus_fumigatus -s sample.fa -t Sample.arf
Question: Hi,

I am trying to create a HISAT2 index with annotation for Rattus_norvegicus (RAT) genome I downloaded from the Ensembl release 94.

I am currently using 220GB memory with 16 cores. My assumption is the memory which I am providing is adequate enough. But I can not create the HISAT2 index, and it gives the error of

Reserving space for joined string
Joining reference sequences
Time to join reference sequences: 00:00:16
Time to read SNPs and splice sites: 00:00:04
is not reverse-deterministic, so reverse-determinize...
Ran out of memory; automatically trying more memory-economical parameters.
is not reverse-deterministic, so reverse-determinize...

and eventually fail with

Could not find approrpiate bmax/dcv settings for building this index.
Switching to a packed string representation.
Total time for call to driver() for forward index: 08:45:56

HISAT2 website does have rat index but they do not have the annotation.

Iam using the command

hisat2-build -p 16 --exon ${EXON} --ss ${SPLICE} ${FASTA_File} ${BASE_NAME}

to create the index.

Any ideas is greatly appreciated.

Thanks

Answer: I think the answer is to provide more memory for the run. Thank you ATpoint.
Question: I'm trying to use fetch to download a sequence. However, I keep on getting `501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)` even though I have ran `sudo cpan install LWP::Protocol::https`

The terminal output is as follows:

$ efetch -db nuccore -id AF086833 -format fasta
501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)
No do_post output returned from 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AF086833&rettype=fasta&retmode=text&edirect=7.00&tool=edirect&email=kjsdhjasv@MacBook-Air.local'
Result of do_post http request is
$VAR1 = bless( {
'_rc' => 501,
'_content' => 'LWP will support https URLs if the LWP::Protocol::https module
is installed.
',
'_msg' => 'Protocol scheme \'https\' is not supported (LWP::Protocol::https not installed)',
'_request' => bless( {
'_method' => 'POST',
'_content' => 'db=nuccore&id=AF086833&rettype=fasta&retmode=text&edirect=7.00&tool=edirect&email=kjsdhjasv@MacBook-Air.local',
'_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
'_headers' => bless( {
'content-type' => 'application/x-www-form-urlencoded',
'user-agent' => 'libwww-perl/6.26'
}, 'HTTP::Headers' )
}, 'HTTP::Request' ),
'_headers' => bless( {
'client-warning' => 'Internal response',
'content-type' => 'text/plain',
'::std_case' => {
'client-warning' => 'Client-Warning',
'client-date' => 'Client-Date'
},
'client-date' => 'Mon, 02 Oct 2017 15:56:40 GMT'
}, 'HTTP::Headers' )
}, 'HTTP::Response' );

Answer: Hi,

How did you install Entrez direct tool?

This is due to a mistake in the installation of Entrez direct tool. Try this code to install it!

mkdir -p ~/src
curl https://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.zip > ~/src/edirect.zip
unzip -o ~/src/edirect.zip -d ~/src
echo 'export PATH=~/src/edirect:$PATH' >> ~/.bashrc
source ~/.bashrc

Tarek
Question: Hi everyone,

I am trying to find some public RNAseq data from single cells. Do you know any of them? Thanks for the help!

Answer: Here you will find one: http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE59739

It is from this publication:
Usoskin, D. et al. Unbiased classification of sensory neuron types by large-scale single-cell RNA sequencing. Nat Neurosci 18, 145–153 (2015).
Question: Hi,

We're working on a RNA-Seq data set, where we used spike In to control for transcriptional amplification. After running the analyses in the standard way and not getting too many significantly regulated genes (~70), we adjusted the normalization according to the spike Ins and got over 2000 genes with significant adj. p-value.

We first calculated the size factors for the Spike In subset and gave it to the DESeq object before running the analysis.

I have a question about the number of reads we are getting after the ercc normalization.

when calculating the sum of reads before and after the ercc normalization we are getting the following numbers

```
sample before after
ctrl1 24762847.36 33353900.56
ctrl2 24973604.51 30803080.63
ctrl3 24727427.27 22561350.5
treat1 24875474.53 18880174.62
treat2 24911409.33 23560508.28
treat3 24588778 21308227.81
```

as you can see there is a clear trend of less reads in the treated samples after taking the spike Ins into account and correcting for the transcriptional bias. This goes perfectly in accord with our expectations.

Our assumption is that the reason for that is that there are less polyA transcripts in the treated data set (an assumption which we would like to verify).

therefore I have two questions -

Are we correct in assuming that the changes in the number of reads is might to less polyA transcripts in the treated samples?

Secondly, is there a statistical way to calculate this diffenrence/change in numbers and give it some kind of a statistical value to it ,so that we can say how significant these results are?

thanks

Assa

Answer: Changes in polyA levels (whether that's due to transcription or turn-over is a different question (enjoy doing the pulse-chase experiment)) are certainly one possibility. There are, of course, others. If this is RNAseq from bulk tissue, then you have to worry about things like ensuring you had highly comparable amounts of starting tissue (in terms of weight, cell number, cell volume, etc.) and that the extraction efficiency is comparable between conditions. If these are single-cell experiments (you should probably have more replicates and have used/counted unique molecular identifiers!), you still have issues of cell size/shape. After all, talking about transcriptional amplification is less interesting if it's just because the cells are huge (or tiny, for the opposite).

Regarding a statistical test, you could just do a T-test between the summed normalized counts. That's less telling (and not significant, given the N of 3) than mentioning ~2000 DE genes, but it's a single number at least. I would suggest looking at the balance of the direction of change of the DE genes. If you really do have a difference due to transcriptional amplification, then the fold-changes should be almost entirely in a single direction.
Question: I want to store a docker image with the input data, code, software, intermediate results and results. This will be a massive file.

Is there something like the sequence archive for people who want to upload a docker image of their whole study? If not, what is the best/cheapest way to store such a file?

Any related thoughts not directly answering the question are also welcome.

Answer: [figshare][1] might be an option. You may need to split the data into chunks and provide a download script, as in [this example][2] (6GB of compressed data split into chunks), though their individual accession size limits may have increased since then.

[1]: https://figshare.com
[2]: https://figshare.com/articles/metaseq_supplemental_data/1091467
Question: Hi friends,

I am going to 8 bacterial bulk RNAseq samples. I am interested in doing following bioinformatics analyses for the 12 samples.

a) Trimming of the reads

b) Alignment of the RNA-seq reads to the reference bacterial genome

c) ID and construction of splice-junctions

d) Reports of known transcripts with annotation and abundance

e) Report of novel transcripts and abundance

f) Testing differential expression

g) GO and Kegg annotation and enrichment analysis

h) Also, DeNovo Assembly

**1) For mammalian, I used the following softwares Tophat, cufflinks, cuffmerge, cuffquant, and cuffdiff?
2) Is there any separate tools/pipelines for bacterial RNAseq data?
3) What are tools for GO, Kegg and enrichment analysis?
4) Which denovo assembly tool is recommended for bacteria?**

Answer: For # 1 since you would not expect splicing to be an issue you could use any aligner (I recommend BBMap suite for trimming/alignment).
#2 there are dedicated pipelines available for bacterial RNAseq. [Rockhopper][1] is a web based one. Search Entrez for others.
#3 You should be able to use GeneSCF (https://www.biostars.org/p/108669/ ).
#4 [SPAdes][2] is about the best bacterial genome assembler.

[1]: https://cs.wellesley.edu/~btjaden/Rockhopper/
[2]: http://bioinf.spbau.ru/en/spades
Question: Hello,

How to define the open chromatin region in ATACseq? The tutorial (https://rockefelleruniversity.github.io/RU_ATAC_Workshop.html) saids insert size less than 100 bp is open chromatin region and insert size in (180, 240) is mono-nucleosome region, I am wondering what is the insert size in (100, 180) belong to?
I know the 9th column of .bam file is the inset size, how to filter the .bam file based on insert size?
I tried
```
samtools view $(FILE)"_mm10.bam" | awk '$9<100' | awk '$9>0' | samtools view -b > openregionbam/$(FILE)"_mm10.bam"
```

But get error
```
[E::sam_parse1] missing SAM header
[W::sam_read1] parse error at line 1
[main_samview] truncated file.
```

What should I do?
Thanks a lot.

Answer: Try this:

It adds `-h` to the first `samtools` command to preserve the header (absence of header is the cause of your error), it preserves that header in the `awk` command (skipping lines starting with `@`) and it respects that alignments to the bottom DNA strand would have negative insert size so you have to filter for insert size being larger than 100 and smaller than -100.

LEN=100
samtools view -h in.bam | awk -v LEN=$LEN '{if ($9 <= LEN && $9 >= -(LEN) && $9 != 0 || $1 ~ /^@/) print $0}' | samtools view -o out.bam

Yes, 100bp is a cutoff I've seen people (and myself) using before but maybe there are (depending on your application) specialized tools that call NFRs mor reliably. NucleoATAC was one that can do it back afair it was quite a pita to use.

Pro-tip: You can replace `awk` by `mawk` to speed up the process a lot.
Question: I can find this description online: QNAME: Query template NAME. Reads/segments having identical QNAME are regarded to come from the same template. A QNAME ‘*’ indicates the information is unavailable. In a SAM file, a read may occupy multiple alignment lines, when its alignment is chimeric or when multiple mappings are given.

The qname example is: M04402:7:000000000-ANGUV:1:2103:23029:8212

This is from a Miseq run, so what is the meaning of these number?

Thanks!!

Answer: If you are only interested in the meaning of that header (referred to as `qname` in SAM format) then check [this link][1].

[1]: https://en.wikipedia.org/wiki/FASTQ_format#Illumina_sequence_identifiers