cannin/biostars_qa huggingface dataset:
Question: I know the Biology behind `deferentially expressed transcripts`.

But I do not understand **what does it mean / difference between** `deferentially expressed genes` and `deferentially expressed transcripts` / isoform levels?

In my RNA-Seq analysis, there are sig. difference data at DEG levels, but those genes are not present at the transcript level.
Is this possible? or my analysis went wrong?

Hypothetical e.g.:

@ At1g16610 ; I have sig results

But @ At1g16610.1 /At1g16610.2 / At1g16610.3 etc nothing is significant.

** I am a wet-lab person, learning bioinformatics thanks to this great community at BIOSTAR

Answer: There are 3 ways changes in expression can occure in a gene with multiple transcripts:

1. All expressed transcripts can be regulated in the same direction (up or down) giving rise to a differential gene expression. In this case the individual transcripts do not necessarily have to be differentially expressed by themself for the gene to be differentially expressed (since the many smaller individual transcript changes when added together becomes a larger change) - but naturally they can be.
2. A subset of expressed transcripts can be regulated. An easy example is that a single transcript is upregulated while the rest are unchanged aka a differentially expressed transcript. This could also lead to a differentially expressed gene if the transcript is one of the most expressed in that gene or the change is very large. But if the transcript is lowly expressed (compared to the other transcripts) this might not cause the gene to be differentially expressed..
3. The third possibility is an isoform switch where one transcript is upregulated while another transcript is downregulated. In this case the transcripts might be differentially expressed - but since they cancel each other out on a gene level the gene can remain unchanged.

It is however important to remember that these events are not mutually exclusive. You can easily have an isoform switch in a gene which is also differentially expressed.

Hope this helped.

Cheers
Kristoffer

P.s. if you are interested in isoform switches my R package [IsoformSwitchAnalyzeR][1] is made to identify isoform switches with predicted functional consequences both in individual genes and on a genome wide level - and it works directly from Salmon/Kallisto data as well. Via these links you can see examples of the [gene level analysis][2] and [genome wide level analysis][3].

[1]: https://bioconductor.org/packages/devel/bioc/html/IsoformSwitchAnalyzeR.html
[2]: https://bioconductor.org/packages/devel/bioc/vignettes/IsoformSwitchAnalyzeR/inst/doc/IsoformSwitchAnalyzeR.html#the-switchplot
[3]: https://bioconductor.org/packages/devel/bioc/vignettes/IsoformSwitchAnalyzeR/inst/doc/IsoformSwitchAnalyzeR.html#genome-wide-summaries
Question: Hi,

I'm interested in to obtain the **external gene names** and the respective **description** based on gene ids for the fruit fly (D. melanogaster) model. I thought that `biomaRt` R package would be a good solution. Therefore I tried the following:

## Convert Geneid into GeneName
#
ensembl <- biomaRt::useMart(biomart = "ensembl",
dataset = "dmelanogaster_gene_ensembl") # select ref db: Drosophila melanogaster genes (BDGP6.28) BDGP6.28
geneName <- biomaRt::getBM( attributes = c("ensembl_gene_id", "external_gene_name", "description"),
filters = "ensembl_gene_id",
values = data$Geneid,
mart = ensembl )

Where the input list is (`data$Geneid %>% head`):

[1] "FBgn0031081" "FBgn0052826" "FBgn0031085" "FBgn0062565" "FBgn0031088" "FBgn0041626"

And the result is (`head(geneName)`):

ensembl_gene_id external_gene_name description
1 FBgn0000022 ac NA
2 FBgn0000045 Act79B NA
3 FBgn0000061 al NA
4 FBgn0000092 CG2650 NA
5 FBgn0000100 RpLP0 NA
6 FBgn0000108 Appl NA

The **external gene name** is retrieved, but not the **description**.

Anyone knows if it is possible to get the **description** from *D. melanogaster* ensembl gene ids with the `biomaRt` R package?

If not, can you recommend a R package tool or any other programmatic tool (rather than go to the site and submit a huge list of genes)?

Thanks in advance,

António

Answer: Using [Entrezdirect][1]:

$ more test
FBgn0000022
FBgn0000045
FBgn0000061

$ for i in `cat ./test`; do printf ${i}"\t"; esearch -db gene -query ${i} | esummary | xtract -pattern DocumentSummary -element NomenclatureSymbol,NomenclatureName; done
FBgn0000022	ac	achaete
FBgn0000045	Act79B	Actin 79B
FBgn0000061	al	aristaless

[1]: http://bit.ly/entrez-direct
Question: I tried to extract reads-pairs aligned concordantly exactly 1 time from Bowtie2 (version 2.3.4.1 or high).
This was fixed in Bowtie2's version a few years ago ([the previous solution][1] and [another solution][2]), but now the problem has resurfaced.

When I run:

bowtie2 -x sample.index -p 50 -t -1 sample_R1.trimmed.fastq -2 sample_R2.trimmed.fastq -S example.sam -k2 --no-discordant --no-mixed

The results would look like this:

Time loading reference: 00:00:00
Time loading forward index: 00:00:00
Time loading mirror index: 00:00:00
Multiseed full-index search: 00:00:35
20826672 reads; of these:
20826672 (100.00%) were paired; of these:
15275620 (73.35%) aligned concordantly 0 times
5048169 (24.24%) aligned concordantly exactly 1 time
502883 (2.41%) aligned concordantly >1 times
26.65% overall alignment rate
Time searching: 00:00:35
Overall time: 00:00:35

What I want to separate and extract is the '5048169 (24.24%) aligned concordantly exactly 1 time' reads-pairs.

But when I followed the previous workarounds, I found:

*1)*

samtools view -hf 0x2 example.sam | grep -v "XS:i:" | wc -l
or

samtools view -bS example.sam > example.bam; samtools view -hf 0x2 example.bam | grep -v "XS:i:" | wc -l

and result(s): `10196512`, which is not equal to 5,048,169*2 = 10,096,338 reads that I want.

*2)*

samtools view -hf 0x2 example.sam | grep -v "XS:i:" | foo.py | wc -l
or

samtools view -bS example.sam > example.bam; samtools view -hf 0x2 example.bam | grep -v "XS:i:" | foo.py | wc -l

and result(s): `10193259`, which is not equal to 5,048,169*2 = 10,096,338 reads that I want.

`foo.py` is (from [Mr. Devon Ryan's solution][1]):

#!/usr/bin/env python
import csv
import sys

f = csv.reader(sys.stdin, dialect="excel-tab")
of = csv.writer(sys.stdout, dialect="excel-tab")
last_read = None
for line in f :
#take care of the header
if(line[0][0] == "@") :
of.writerow(line)
continue

if(last_read == None) :
last_read = line
else :
if(last_read[0] == line[0]) :
of.writerow(last_read)
of.writerow(line)
last_read = None
else :
last_read = line

*3)*

samtools view example.sam | grep YT:Z:CP | uniq -u | wc -l
or

samtools view -bS example.sam | grep YT:Z:CP | uniq -u | wc -l

and result(s): `12107870`, which is not equal to 5,048,169*2 = 10,096,338 reads that I want.

*4)*

samtools view -hf 0x2 example.sam | grep -v "XS:i:" | grep YT:Z:CP | uniq -u | wc -l
or
samtools view -hf 0x2 example.sam | grep -v "XS:i:" | grep YT:Z:CP | wc -l
or
samtools view -hf 0x2 example.bam | grep -v "XS:i:" | grep YT:Z:CP | wc -l

and result(s): `10159819`, which is not equal to 5,048,169*2 = 10,096,338 reads that I want.

**I've tried almost every existing solution but to no avail. Urgently ask for your help.**

[1]: https://www.biostars.org/p/95929/#9536541
[2]: https://github.com/BenLangmead/bowtie2/issues/58

Answer: This might help.
https://github.com/BenLangmead/bowtie2/issues/406
Question: When I run the script below, output is getting split into single chars. Any idea why? It looks like the second argument gets split into single chars. I am trying to align the word sequences. I will have many words hence cannot map them to letters only.

Thx Ercan

from Bio.Seq import Seq

from Bio.pairwise2 import format_alignment

fruits = ["orange","pear", "apple","pear","orange"]

fruits1 = ["pear","apple"]

from Bio import pairwise2

alignments = pairwise2.align.localms(fruits,fruits1,2,-1,-0.5,-0.1, gap_char=["-"])

for a in alignments:

print(format_alignment(*a))

Output:

['orange', 'r', 'a', 'e', 'p', 'e', 'l', 'p', 'p', 'a', 'pear', 'orange']
|||||||||
['-', 'r', 'a', 'e', 'p', 'e', 'l', 'p', 'p', 'a', '-', '-']
Score=4

Answer: BioPython allows for performing an alignment on sequences' symbols (not words). I can't think of any program that would do the words' alignment. I suggest you search the Web for Python modules/functions that perform string character alignments and then modify the existing code to handle words.

For example, I tweaked [some code from GitHub][1] to do local (*water*) and global (*needle*) alignments of words.

Filename: `alignment.py`:

https://gist.github.com/aziele/6192a38862ce569fe1b9cbe377339fbe

**Usage:**

import alignment
fruits1 = ["orange", "pear", "apple", "pear", "orange"]
fruits2 = ["pear", "apple"]
aln = alignment.needle(fruits1, fruits2)
identity = aln[0]
score = aln[1]
print(identity, score)
print(aln[2])
print(aln[3])

Output:

(0.4, 25)
['-', 'pear', 'apple', '-', '-']
['orange', 'pear', 'apple', 'pear', 'orange']

[1]: https://github.com/alevchuk/pairwise-alignment-in-python/blob/master/alignment.py
Question: Hi,
I would like to know the AF of ~4000 rsIDs (dbSNP) quickly in Gnomad GRCh37 database. I looked at the instructions for getting this from Gnomad and it seems I would need to download the entire data. Is there is a faster and less involved way fo getting the AF from Gnomad without downloading more than 50GB of data? Any script or API or another way to crawl Gnomad database?
Thanks!

Answer: The `af-only-gnomad.raw.vcf` file [here](https://console.cloud.google.com/storage/browser/gatk-best-practices/somatic-b37) is 13 GB. It might work better than the 80+ GB file from the gnomAD website.
Question: I've been dealing with comparison between 40 bed files to one bed file using `intersectBed -a -b` command. I'm just wondering, is there any commands in [Bedtools][1] which can help us compare multiple bed files?

Say, I have 40 bed files and one particular bed file. I want to identify those regions in the 40 bed files overlaps with the particular one. I mean the 40 to 1 comparison.

Is there any fast ways to compare them and do not need type the code one by one like `intersectBed -a -b`, `intersectBed -a -c`, `intersectBed -a -d`, `intersectBed -a -e`...

[1]: https://github.com/arq5x/bedtools2

Answer: [One fast option is to use BEDOPS][1]:

$ bedops --intersect A.bed B.bed C.bed ... > answer.bed

You can use [lots of inputs][2] efficiently. The input files just have to be [sorted][3].

The above command intersects `B.bed C.bed` etc. with `A.bed`, reporting all elements of `A.bed` that overlap `B.bed C.bed` etc.

Let's say you want to go the other direction efficiently. You can use BEDOPS with UNIX pipes and redirect standard output from one command to the next:

$ bedops --everything B.bed C.bed ... | bedops --intersect - A.bed > answer.bed

This does a **multiset union** of all the elements in `B.bed C.bed` etc. and passes these to an `--intersect` operation with `A.bed`.

The result file reports all elements of `B.bed C.bed` etc., which overlap `A.bed`.

The difference between these two directions is in which sets of elements get reported in the overlap. In the first case, elements of `A.bed` are reported. In the second case, elements of `B.bed C.bed` etc. are reported. Generally, this is not a symmetric operation.

If you have a lot of files to sort, a quick bash one-liner can take care of this:

$ for fn in `ls *.bed`; do sort-bed ${fn} > ${fn%.*}.sorted.bed; done

Some use GNU `sort` to do sorting of BED files, but BEDOPS `sort-bed` is usually faster.

[1]: http://bedops.readthedocs.org/en/latest/content/usage-examples/multiple-inputs.html#multiple-inputs
[2]: http://bedops.readthedocs.org/en/latest/content/reference/set-operations/bedops.html#bedops
[3]: http://bedops.readthedocs.org/en/latest/content/reference/file-management/sorting/sort-bed.html
Question: I am working with certain codes in python and I wanted to know if a certain value of a variable can be accesssed outside the function without explicitly declaring it a global variable.
for example:

I was using this:

def func1():
v1 = "value1"
return v1
print ("This is the first function")

def func2():
v2 = "value2"
func1()
print ("This is the second function")

fo = input ("Enter 1/2\n")
if fo == '1':
func1()
elif fo == '2':
func2()
else:
print ("invalid")
print (v1)

It results me this error:

print (v1)
NameError: name 'v1' is not defined

Can I use this value of `'v1'` variable outside the function without defining it to b `'global'`??

Answer: You have two options (actually more but lets keep it with two)

As Devon Ryan suggest and also my preference:

def func1():
v1 = "value1"
print ("This is the first function")
return v1

def func2():
v2 = "value2"
func1()
print ("This is the second function")

fo = input ("Enter 1/2\n")
if fo == '1':
funcValue = func1()
elif fo == '2':
funcValue = func2()
else:
print ("invalid")

print (funcValue)

Also you need to print first before returning something. After the return the function stops.

Or make a global variable, but don't do this if it is not necessary:

v1 = ""
v2 = ""
def func1():
global v1
v1 = "value1"
print ("This is the first function")
return v1

def func2():
global v2
v2 = "value2"
func1()
print ("This is the second function")

fo = input ("Enter 1/2\n")
if fo == '1':
func1()
elif fo == '2':
func2()
else:
print ("invalid")

print (v1)
Question: Where is the best standard place to store single-cell datasets to be accessible for the community?
Is GEO the right place? I looked up some single-cell datasets and the number of samples is quite low does each sample include all the cells fastq files?

- [GSE57249][1]: Has just 56 samples
- [GSE57872][2]: Has 587 samples
- [GSE145926][3]: Has 21 samples
- ?

Related to sharing the data set, is it worth to publish a Data Note of a single cell data set? (Tradeoff between salami-splicing publications and work on a new manuscript, citation/reference of the data set for multiple purposes...)

[1]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE57249
[2]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE57872
[3]: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE145926

Answer: I was running on the internet with the same query and I think yes GEO is the best place to deposit raw sequencing data. There are other single cell specific data repositories, for example, Broad Institute's [Single Cell Portal][1], [UCSC Cell Browser][2], [PanglaoDB][3], [Singel Cell Studies][4], and perhaps plenty of others. Usually, these portal has an advantage over GEO, where you can interactively access the scRNA-Seq data like visualizing cluster embeddings, visualizing gene expression, cluster-specific pathways distribution, etc. Nevertheless, you can upload the processed data (expression matrices, metadata, R objects, etc) to GEO, but the user needs to download and visualize it by themselves.

I have uploaded 10X genomics data to the GEO once and in my opinion, you can directly submit the samples (multiplexed; raw sequencing data) so that users can start processing your data from the scratch.

So in my case, we had 3 samples with approximately 14k cells in each. I submitted 3 samples with the processed data (by CellRanger, feature count matrix, metadata, cell barcodes, and Seurat object in .rds format).

I hope this would help you during your submission.

Cheers!!

;)

[1]: https://singlecell.broadinstitute.org/single_cell
[2]: https://cells.ucsc.edu/
[3]: https://panglaodb.se/papers.html
[4]: https://www.nxn.se/single-cell-studies/
Question: Dear all,

I have WES bam files from tumour samples and a normal sample of one patient, and I would like to estimate the tumour purity of the tumour samples with ABSOLUTE. They say that "you can supply a tab delimited segmentation file (e.g. from array CGH or massively parralel sequencing experiments) - this file must contain the columns "Chromosome", "Start", "End", "Num_Probes" and "Segment_Mean".".

I have read similar questions here and here, but I still haven't understood how to create the segmentation file.

I have already analysed copy number alterations (CNAs) using CopywriteR and it gives me a read_counts file with "Chromosomes", "Start", "End" and the read counts for the samples but I don't see where to get the "Num_probes" and the "Segment_Mean" information.

I would highly appreciate your help

Answer: If you've run all 3 phases of the CopywriteR pipeline (the functions `preCopywriteR`, `CopywriteR`, and `plotCNA`) according to the tutorial/vignette/manual, then in the output folder, under the `CNAprofiles` subdirectory where you find `read_counts.txt` and all the other output files, there will be a file called `segment.Rdata` -- this is generated by `plotCNA`.

If you load `segment.Rdata` in R, you get a list with an `$output` field that contains the dataframe you need, essentially in SEG format. Here's my script for loading that file and converting it to SEG:

```r
#!/usr/bin/env Rscript
# Usage: Rscript copywriter2seg.R

(args = commandArgs(TRUE))

in_fname = args[1] # e.g. segment.Rdata
out_fname = args[2] # e.g. sample.seg
cat("Loading", in_fname, "\n")
load(in_fname)
seg = segment.CNA.object$output
seg$loc.start = floor(seg$loc.start) # force coordinates to be integers
seg$loc.end = floor(seg$loc.end)
write.table(seg, out_fname, quote=FALSE, sep="\t", row.names=FALSE, col.names=TRUE)
cat("Wrote", out_fname, "\n")
```
Question: When I try this [vcftools_results]$ /home/X/vcftools_0.1.6/perl/vcf-validator sample1.flt.vcf.gz

it gives an error.

Can't locate Vcf.pm in @INC (@INC contains: /share/lib64/perl5/5.8.8 /share/lib64/perl5 /share/lib/perl5/5.10.0/x86_64-linux-thread-multi-ld /share/lib/perl5/5.10.0

/share/lib/perl5/5.8.8 /share/lib/perl5 /usr/lib/perl5/5.8.8 /usr/lib/perl5

/share/lib/perl5/5.10.0/x86_64-linux-thread-multi-ld /share/lib/perl5/5.10.0

/share/lib/perl5/site_perl/5.10.0/x86_64-linux-thread-multi-ld

/share/lib/perl5/site_perl/5.10.0 /share/lib/perl5/site_perl/5.8.8

/share/lib/perl5/site_perl .) at /home/X/vcftools_0.1.6/perl/vcf-validator line 9.

BEGIN failed--compilation aborted at /home/X/vcftools_0.1.6/perl/vcf-validator line 9.

how can I fix it?

Answer: I'm guessing that Vcf.pm is located somewhere under /home/X/vcftools_0.1.6/perl.

So the problem is that this is not a standard Perl library location (@INC) and Perl cannot find it.

Did you read the documentation? It contains this important piece of information:

For running the Perl scripts, the PERL5LIB environment variable must be set to include the Vcf.pm module

export PERL5LIB=/path/to/your/installation/perl
Question: Hi,

I know that there have been multiple suggestions for getting gene ids. I have used a method to do it myself and this is just a post to see if I am doing it right. I have the oligonucleotide probe ids for a set of transcription factors and genes which are part of the Apis Mellifera gene regulatory network. What I did was:

1) I went on GEO and retrieved the sequences for the probes.

2) I further went on NCBI and used BLAST to match these probe sequences with all the sequences in the database related to *Apis Mellifera*

3) Further, I extracted the ENTREZ IDs of all the hits in an excel sheet.

4) I converted the ENTREZ IDs for all the hits into UNIPROT IDs on UNIPROT.

5) I further, used these UNIPROT IDs to do a functional annotation analysis on DAVID.

Please advice if the methodology I followed was correct and if not, please suggest any changes that I should make to it.

Answer: So, basically, since there could be a number of hits for a sequence potentially leading to unwanted results, I found a list mapping the probes with their ENTREZ ids and used that for Analysis on DAVID. If you have no option other than using sequences of probes, I would suggest using MEGABLAST to minimize matches. This worked for me and there were only minimal matches per sequence, but I figured that using the ENTREZ ids would be a better option for me. Another suggestion which I received and would like to share is to be careful about using Excel. Because Excel could potentially change ids. Hope these suggestions help you with your research.
Question: Hi everyone,

I'm trying to find the genomic location of a list of sequences in the reference genome.
So if my sequence is 'ACGTACGTAGTCATGC', I want as an output something like this:

chr1 position of the first nucleotide position of the last nucleotide

Are there any tools to do this? Blastn maybe? But since I have the reference genome locally installed, there might be some other options?

Answer: If you only need exact matches, then give one of the answers in this thread a try: https://www.biostars.org/p/66752/
Question: I have 12 pair-end samples, I would like to compare just Fpkm values using tophat/cufflinks. Would the Fpkm values be the same (for annotated and/or novel transcripts) regardless whether or not I process all 12 samples in the same run (i.e., put 12 group folders into the main folder in which I run the program) or if I run 2 samples at a time and then just compile all 12 in one table?

Could I process and compare in the same run pair-end samples sequenced 50bp or 100bp from each end? Would the platform affect the Fpkm values or I can assume normalization is appropriate and I can rely on differential expression results?

Answer: They won't be the same if you have novel transcripts and then use multiple runs. There's a step where one merges (via cuffmerge) the annotations produced by cufflinks for each sample. The FPKMs are dependent upon that merged annotation, so at least if you perform the merge step with different groups of samples then you're pretty much guaranteed to see differences. Further, I think cuffdiff defaults to using a DESeq-like library size normalization rather than the old and less reliable "divide by the number of aligned reads" method. That can also produce (generally pretty small) differences in FPKM estimation.

Regarding using a mix of read lengths, yeah, you can do that (you can also mix paired-end and single-end). Keep in mind that mappability will change a bit, so that'll produce a bias. If you have multiple platforms, then that'll be a confounder. I should note that cuffdiff is really only designed to handle very simple experimental designs (i.e., two group comparisons or a time series). If your experiment doesn't fit into one of the defaults, you'd be better off using a different tool (just use one of the count-based methods, like DESeq2/edgeR/etc.).
Question: Hello,

I am trying to find a qc package to be used with.bam files. In using the RNA-SeQC package I ran into difficulty because it required -s, -r, and -r arguments. I only have a .bam file and do not have these required parameters.

Is there a RNA-Seq specific QC tool where I can simply point to the .bam file and have it create a qc report?

Thank you.

Answer: [RSeQC](http://rseqc.sourceforge.net/), [QoRTs](https://hartleys.github.io/QoRTs/), [BamQC](https://github.com/s-andrews/BamQC) are all valid options.

However, I'm surprised you cannot use RNA-seQC, because in its description it says: "RNA-SeQC is a java program which computes a series of quality control metrics for RNA-seq data. The input can be one or more BAM files."