cannin/biostars_qa huggingface dataset:
Question: I would like to extract the raw counts from my RNA-Seq experiment from Cufflinks. I know the arguments about alternate splicing but I am trying to compare different differential expression methods in my model plant.

While the Cufflinks manual is very clear. In the Cuffnorm section it says these reads "should not be used with downstream differential expression tools that require **raw** counts as input".

However, the CummeRbund release schedule says that "raw and normalized count tables and associated statistics all features" have been available for some time. From what I can tell this is the obtained from the `count()` function. However in the example the count column values seem normalised. So I would like to check are these values appropriate to use in applications like DESeq and edgeR?

If not is there another way to get these values?

Answer: Just use htseq-count or featureCounts. If some random value column from cufflinks agrees with these values then you'll know for the future which one to use. Cufflinks is a complicated enough program that one needs to go through the code to actually know what's truly going on, so I would strongly dissuade you from using any of it's values in DESeq2/edgeR/etc. without checking in a couple samples if they match what's expected.
Question: Hi all,

Im using Hisat2 to align RNAseq data and encountering following error:

(ERR): hisat2-align died with signal 8 (FPE) (core dumped)

using following command

./hisat2 --mp 0,0 --np 0 --pen-noncansplice 10000000 --min-intronlen 20 --max-intronlen 10000 --novel-splicesite-outfile /data/memona/hisat2-2.1.0/result/ --rna-strandness RF --dta -p64 --summary-file -x /data/memona/hisat2-2.1.0/hisat_index -1 /data/memona/Trimmomatic-0.36/SRR9591_E_1P.fq -2 /data/memona/Trimmomatic-0.36/SRR9591_E_2P.fq -S /data/memona/hisat2-2.1.0/result/hisat_align.sam

can anyone kindly suggest possible cause and solution to the problem.

Thanks

Answer: Looks like these issues have been already taken care in the recent version. Check your HISAT version using

hisat2 --version

Update if necessary
Question: Hello everyone,

It is possible to visualize the progress of a local blast+ (verbose mode or something [maybe python ou perl])?

Thanks in advance!

Answer: I don't think so, but if I'm wrong, I would be happy because I miss it too. What I'm doing now is (works with a tab delimited output and a multifasta file query) to find the last query ID written in the blast output, and then to find it into the fasta file. This works only because sequences are "blasted" in the same order as in the query file.

```
## 1st step
tail -1 out.tab.blast | cut -f 1 > lastqueryblasted

## 2nd step
grep -n -f lastqueryblasted query.fasta
# res1, ex: 23000

## last step, pourcentage computing
wc -l query.fasta
# res2, ex: 33000

# your (under)estimated pourcentage
echo "(23000 / 33000) *100" | bc -l
```

This is not very convenient I know, but it helps.

EDIT: I use that so much often, that I eventually wrapped it in a script.

If you name it `blast_monitor.sh`, it works like this:

blastab_monitor.sh blastresult.blast queryfasta.fa

Here it is:

```
#!/bin/sh
### monitor a tab-outputed blast job by giving the approximative % done
blast=$1
query=$2
echo "the blast out is: "$blast
echo "the fasta query is: "$query
echo
curquery=$(tail -1 $blast | cut -f 1)
curline=$(fgrep -n $curquery $query | cut -f 1 -d ':')
nblines=$(wc -l $query | cut -f 1 -d " ")
percent=$(echo "($curline/$nblines) *100" | bc -l | cut -c 1-4)
echo "The blast job is about $percent % done..."
```

I hope that helps.
Question: Hello,
I have a fasta file with lines of format -

>FBti0019256 type=transposable_element; loc=2L:22300300..22304444; name=invader2{}555; dbxref=FlyBase_Annotation_IDs:TE19256,FlyBase:FBti0019256; MD5=d9259a0e33aad699215e64916bd47a5b; length=4145; release=r6.19; species=Dmel;

I would like to convert these lines into a bed file of format -

chr2L /t 22300300 /t 22304444 /t invader2

Is there a program that can directly perform this conversion or an awk command that can do this easier? Please let me know, thank you for your help.

Answer: My approach is lazier... Just split on the punctuation, pull out the columns you want, and separate by "\t".

awk -F'[ *:.]' -v OFS="\t" '{print "chr"$5,$6,$8,$13}' tmp.txt

# Output
chr2L	22300300	22304444	invader2
Question: Hi!

I am trying to make smaller many BAM files (around 60 of them) of size ~200GB (due to disk space limitations) by removing base qualities and tags and other unwanted information. Doing copy number analysis, for me base qualities and tags and duplicates are somehow unwanted information. What I only care about is the mapping quality (MAPQ) since I filter low quality reads!

Currently, I am using [bamUtils squeeze][1] command. I don't know yet how good this tool is in making the bam file smaller! The *squeeze* sub command can replace QNAME with an integer, remove duplicates, and remove OQ tag (i.e. original base qualities) but not the QUAL field. However, for the QUAL field, the tool provides the binning option (to reduce the number of possible quality scores).

Previously, I used [cgat bam2bam method=strip-quality][2] which deletes only the QUAL field. This tool is slow (takes ~12 hours for a 160 GB BAM file) and didn't free much space. The modified BAM file was only 3GB smaller for a 160GB file.

I was wondering if deleting whatever comes after the SEQ (in a SAM/BAM file) will work (i.e QUAL and all other tags)? and if yes, what would be the fastest way to apply that? Or, if there is a tool available that I was not able to find?

Thanks in advance for sharing your ideas!

EDIT1: My question might have been misleading since I said "I only care about MAPQ". I also care about the FLAG and SEQ. Since later in the pipeline, I will call variants; but there only FLAG, MAPQ and SEQ are needed and not any thing else.

EDIT2: Now, I have the result from using [bamUtils squeeze][1] and to me the result is satisfactory. The BAM file is ~4-fold smaller when one:

1. removes the OQ tag,
2. removes the duplicates, and
3. bin the base quality scores

And it took less than 4 hours (3:51) to squeeze a 160GB file to 43.5GB.

[1]: http://genome.sph.umich.edu/wiki/BamUtil:_squeeze
[2]: https://www.cgat.org/downloads/public/cgat/documentation/scripts/bam2bam.html

Answer: I quickly wrote a tool: see https://github.com/lindenb/jvarkit/wiki/Biostar173114

I checked, it's compatible with samtools.

$ java -jar dist/biostar173114.jar --formatout sam my.bam | samtools view -h | head -n 100

@HD VN:1.5 GO:none SO:coordinate
@SQ SN:1 LN:249250621
@SQ SN:2 LN:243199373
(...)
@SQ SN:GL000192.1 LN:547496
R0 0 1 822010 55 51M * 0 0 * *
R1 16 1 822010 55 51M * 0 0 * *
R2 0 1 829084 60 87M13S * 0 0 * *
R3 0 1 829084 60 87M13S * 0 0 * *
R4 0 1 829084 60 87M13S * 0 0 * *
R5 16 1 829122 59 48M36D52M * 0 0 * *
R6 16 1 829122 59 48M36D52M * 0 0 * *
R7 16 1 829122 59 48M36D52M * 0 0 * *
R8 0 1 848195 46 19S20M3S * 0 0 * *
R9 16 1 848195 46 19S20M3S * 0 0 * *
Ra 0 1 908830 60 47M * 0 0 * *
Rb 16 1 908830 60 47M * 0 0 * *
Rc 0 1 922930 60 100M * 0 0 * *
Rd 0 1 922930 60 100M * 0 0 * *
Re 0 1 922930 60 100M * 0 0 * *

and for better results, you should process your data in parallel.
Question: I have a vcf file with contains GT:DS:GP for FORMAT field, how can i only keep the GT?
The vcftools "--extract-FORMAT-info GT" command will at the same time filtering all other columns. But i wish to keep those.

Answer: > https://samtools.github.io/bcftools/bcftools.html
>
>... and "FORMAT" to remove all FORMAT tags except GT.

bcftools annotate -x 'FORMAT' in.vcf
Question: I want the total lengths of regions that are exons, splicing, UTRs, and ncRNA. I can't distinguish them from the provided bed file. If there is a way using the bed file, please let me know. Thank you.

Answer: You can use a GTF with annotated intervals for the categories you want. For example this [one][1] from Ensembl. Use something like unix grep to create files for the categories you are interested. Then use [bedtools intersect][2] to get the regions of your BED file overlapping with each category GTF file. Finally use the procedure described [here][3] to find total lengths.

Also, if you are only interested on the ratio for each category, you could simplify this by using [bedtools jaccard][4]. In this case you would only need your BED file and the categories GTF files.

[1]: ftp://ftp.ensembl.org/pub/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh38.84.gtf.gz
[2]: http://bedtools.readthedocs.io/en/latest/content/tools/intersect.html
[3]: https://www.biostars.org/p/68283/
[4]: http://bedtools.readthedocs.io/en/latest/content/tools/jaccard.html
Question: I am just starting to learn to use bioinformatics tools. My university has a limited and expensive bioinformatics team, so I'm mostly on my own except for big questions.

I am planning to use GATK to run 58 cancer control/normal pairs of Exome sequencing data (Illumina) from FASTQ or BAM file format, through the pipeline, with an output of a VCF & MAF format for analysis.

The current GATK pipeline is used for disease but not cancer, so I was wondering if anyone knew if there should be changes made for cancer. Here's the current pipeline starting with BAM files:

(Non-GATK) Picard Mark Duplicates or Samtools roundup

Indel Realignment (Realigner TargetCreator + Indel Realigner)

Base Quality Score Reacalibration (Base Recalibrator + PrintReads)

HaplotypeCaller <- I've been told this is for germline variants; what can I use for somatic variants?

VQSR (VariantRecalibrator and ApplyRecalibrator in SNP and INDEL mode)

Annotation using Oncotator (?)

I'd like some verification that this pipeline will output what I need to run my samples on MuTect, MutSig, or some other analysis program. I appreciate any advice.

Crossposted on Stack Exchange Biology.

Answer: We (GATK docs team) are working on some docs for the somatic variant calling use case. In a nutshell, you'll need to do an additional pre-processing step called co-cleaning where you perform indel realignment on the tumor and normal in a pair together, use ContEst to estimate cross-sample contamination, use MuTect to call variants (not HC, which is not able to call low-AF variants like MuTect), do some manual filtering and processing to eliminate artifacts (VQSR is not appropriate for somatic calls) and finally annotate with Oncotator.

We're happy to provide more details in the GATK support forum. The CGA homepage (http://www.broadinstitute.org/cancer/cga/Home) is also a good resource.
Question: I am mapping a series of fastq files, containing ChiP-Seq data, with bowtie2 (my mapping software of choice) to the GRCh37 genome: `https://www.gencodegenes.org/human/release_39lift37.html` . I have mapped as follows:

bwa mem GRCh37.primary_assembly.genome.fa ${srr}.fastq.gz | samtools view -S -b -q 10 - > ${srr}.bam

I then create the fasta file for the bt2 reference and index it:

samtools faidx GRCh37.primary_assembly.genome.fa
java -jar ~/picard.jar CreateSequenceDictionary -R GRCh37.primary_assembly.genome.fa

However, when I then try to call variants in this data I get an error:

java -jar ~/picard.jar MarkDuplicates -I ${srr}_sorted.bam -M metrics.txt -O ${srr}_sorted_dedup.bam
gatk --java-options "-Xms10G -Xmx10G -XX:ParallelGCThreads=2" BaseRecalibrator \
-I ${srr}_sorted_dedup.bam \
-R GRCh37.primary_assembly.genome.fa \
-O ${srr}_sorted_dedup.report \
--known-sites variants.vcf

The marking of the duplicates part seems to work fine but for the recalibrating the bases part the contigs do not seem to match up:

A USER ERROR has occurred: Input files reference and features have incompatible contigs: Found

contigs with the same name but different lengths:
contig reference = 10 / 135534747
contig features = 10 / 135284542.
reference contigs = [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1, 20, 21, 22, 2, 3, 4, 5, 6, 7, 8, 9, MT, X, Y]
features contigs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, X, Y, XY, M]

I have read that the vcf file that I am using was made according to GRCh37: https://www.synapse.org/#!Synapse:syn10901595

I downloaded the plink files from https://www.synapse.org/#!Synapse:syn17008939 and made the vcf file as follows:

plink --bfile ROSMAP_data --recode vcf --output-chr M

Hence I am wondering how the vcf file can mismatch the ref genome when they are the same- has anyone experienced this problem before?

Also just check that the ref genome should not be b37 does anyone know of a link to the b37 fasta file? The following link which I tried- leads to a NO PERMISSION page for me:

https://console.cloud.google.com/storage/browser/_details/broad-references/hg19/v0/Homo_sapiens_assembly19.fasta;tab=live_object

Answer: I solved this problem by finding the correct liftover chain file to convert my `vcf/known-sites` file to a version that was compatible with my reference genome. The README that described the making of this vcf file states that it was created using the GRCh37 genome as the reference- when in fact it was the b37 genome- hence:

java -jar picard.jar LiftoverVcf \\
I=variants.vcf \\
O=lifted_over.vcf \\
CHAIN=b37tohg19.chain \\
REJECT=rejected_variants.vcf \\
R=GRCh37.primary_assembly.genome.fa

And then the base recalibration step worked.

NOTE: The chain file is downloaded from https://github.com/broadgsa/gatk/tree/master/public/chainFiles
Question: Hi everybody,

I am running GATK pipeline for variant calling in trinity with following command.

/usr/local/bin/Trinityrnaseq-v2.6.6/Analysis/SuperTranscripts/AllelicVariants/run_variant_calling.py --st_fa ./SuperDuper.fasta --st_gtf ./SuperDuper.gff -p ./FCHG1.fq.gz ./FCHG.fq.gz -o ./variant_calls_outdir

But I got following error:

Error, missing path to Picard-Tools in $PICARD_HOME.
based on the instructions on the trinity web, I need to set the installation directory to the environmental variable ${PICARD_HOME}. for this purpose (adding the picard path to $PATH) I tried following command but my problem was not solved.

PATH="$PATH:$HOME/picard-2.23.0"
echo $PATH
output

/home/ubuntu/bin:/home/ubuntu/.local/bin:/home/ubuntu/anaconda3/bin:/home/ubuntu/perl5/bin:/usr/local/sbin:/usr/local/bin:/usr/sbin:/usr/bin:/sbin:/bin:/usr/games:/usr/local/games:/snap/bin:/home/ubuntu/picard-2.23.0

I don’t know how can I set the installation directory to the environmental variable ${PICARD_HOME}.

best regards

Answer: The script is asking for the $PICARD_HOME variable, then you need to create it:

export PICARD_HOME="$HOME/picard-2.23.0"
Question: Hi!

I am currently working to analyze alpha & beta diversity with microbiome data of 48 sample.
(Animal stool samples)

I finished working of Beta diversity and got the data of significance. (Weighted UniFrac, Bray-curtis)

However, I received this message and need some help from you.

The message is as below (From reviewer)
------
I highly recommend the use of **linear modeling (LM) or generalized linear modeling (GLM) which is** commonly used in microbiome studies **rather than a Wilcoxon rank-sum test**. This will allow you to better control for things that may impact your results such as the age animals. It is also important to control for site as this is known to significantly affect the gut microbiome and could be included as a random variable in a generalized linear mixed model (GLMM).
----------

I found that non-binomial regression or possion regression is the most commonly used glm for
microbial analysis and trying to use the models in comapring Alpha and Beta diversity difference of the 16s rRNA data.

However, my question is that what package in R or python is recommended when figuring out microbial
analysis. (Perhaps, many researchers in the metabarcoding field use some common statistical tools)

And, also I want to know that how can I modify (adding a dependent variable or etc. ) the data file
to use in the GLM model. (What variable should be added to compare in Alpha and Beta diversity analysis, and how can I calculate them?)

my metadata file is as below

![Meta data file to be used in GLM][1]

[1]: /media/images/478208ba-66c0-4776-a09d-8228c249

Any recommendation about using Package source, specific GLM model or the way to modify data -
would be a big help for me to solve the problem.

Thank you

Answer: Hello - in addition to the packages listed above (which I have used and they are good), there are also dedicated bioconductor packages written ***specifically for microbiota analysis***.

A few that might work based on your post include: package `DirichletMultinomial`, package `Microbiota`, and package `phyloseq`.

Links to these 3 resources and dozens of other dedicated Bioconductor packages specifically for microbiota [may be found here][1].

[1]: https://www.bioconductor.org/packages/release/BiocViews.html#___Microbiome
Question: Hello,

I have set a UCSC trackhub, which includes annotations that were originally in the GTF or GFF formats. So I followed the steps described here: https://genome.ucsc.edu/goldenPath/help/bigGenePred.html (example 4: GTF (or GFF) to BigGenePred) to obtain a binary bigGenePred for each GTF / GFF.

It works fine, except that when visualizing the result, I loose data on introns and strand orientation. It means that for each gene annotated, I get a solid line, going from the 5' end to the 3'end of the gene.

One example of what I get (the first 3 light blue top lines are from one of my GTF annotation turned to binary bigGenePred):

![Track hub screen capture][1]

I really want to see a finer granularity where you can actually see exons / introns and strand orientation, so I tried to upload one of my GTF to a custom track. This way works: I do not lose data from the annotation and get to see all the details. However, I could not find a way to integrate a custom track to a track hub.

Result with the custom track (top 3 black lines are from one of my GTF annotation):
![Custom track screen capture][2]

So this my question: what is the best way to integrate a GTF file to a trackhub and not lose any detail ?

Thanks for the help.

[1]: https://i.ibb.co/zsqHPKL/Capture-d-e-cran-2020-05-20-a-11-30-05.png
[2]: https://i.ibb.co/mGrVvrD/Capture-d-e-cran-2020-05-20-a-11-45-05.png

Answer: Glad it was solved! I'm just going to copy the response as an answer so the post doesn't show up unanswered:

Hello,

bigGenePred does support intron/strand information. It should not be a problem when converting from a GTF/GFF. Initial suspicions are the problem may be the way the trackDb stanza is declared in the hub. Are you stating type bigGenePred, etc? Another possibility could be the input GTF file or the conversion (stating the correct bidBed 12+8), but that's less likely. Essentially from your first screenshot it looks like the Genome Browser is displaying your file as a BED3.

If you would like to email us a copy of the GTF file as well as a link to the hub to our private mailing list (genome-www@soe.ucsc.edu) we could take a look. Note that only internal Genome Browser staff can see the contents of the message.
Question: Does anybody know how to plot a network similar to the one in the below TCGA paper?

http://www.nature.com/nature/journal/v517/n7536/fig_tab/nature14129_F5.html

More specifically, is there a tool to produce such visualization or was it created manually?

Answer: Maybe you could make it with [OmniGraffle][1]. It has some nice network drawing tools. One thing that OmniGraffle offers that I don't think Illustrator offers is automatic management of connections between elements. For instance, you drag a node in the network and any connected edges are also redrawn. You can align nodes horizontally and vertically, to make a clean layout.

You might also investigate [the SBGN format][2]. One example in the relevant Nature paper shows [stimulation and inhibition][3], for instance. It looks like an [XML-based document format][4] that might lend well to automation, if you are making lots of these figures.

[1]: https://www.omnigroup.com/omnigraffle
[2]: http://www.sbgn.org/Main_Page
[3]: http://www.nature.com/nprot/journal/v7/n3/full/nprot.2012.002.html?WT.ec_id=NPROT-201203
[4]: http://www.sbgn.org/LibSBGN/Exchange_Format
Question: I'm using the EMBOSS `transeq` tool to translate the first ORF of 26,000 sequences. Since the tool is pretty slow (takes ~1 minute per short sequence) and to make the process parallel, I split the fasta file into smaller files (down to one sequence per file) and then run `transeq` on each file. However, the sequences are renamed in the following format: `EMBOSS_001_1`.

How could I prevent `transeq` from renaming these sequences?

If nothing else works, I'll create a script that manages the translation of each individual sequences and makes sure to rename it after it has been translated.

Answer: Assuming that your input sequences are in [fasta sequence format][1] then one option which might help is to disable the parsing of the sequence identifier. By default [EMBOSS][2] attempts to identify the appropriate sequence format from the initial lines of the input. If the input is identified as being in a fasta sequence format variant (i.e. 'fasta' format) then identifier parsing will be used to extract information commonly encoded in structured identifiers, such as the database name, accession and entry name. For a list of possible fasta sequence format variants supported by EMBOSS see the ["EMBOSS Users Guide"][3] appendix ["A.1. Supported Sequence Formats"][4]. The additional processing can be disabled by specifying the generic 'pearson' format fro the input, which uses the first non-whitespace token on the header line as the sequence identifier. The input sequence format is specified using the `-sformat` option, so your command-line will look something like:

transeq -frame 1 -sformat pearson -sequence inSeqFile.fa -outseq outSeqFile.fa

Please note that the EMBOSS 6.6.0 release contains a bug which makes the parsing of sequence identifiers consisting of two pipe separated fields (e.g. `db|seqId`) very slow, which might explain the poor performance you are seeing. A fix for this problem should be available in the post release patches (see ftp://emboss.open-bio.org/pub/EMBOSS/fixes/), and the issue did not appear in the previous release (EMBOSS 6.5.7).

[1]: http://en.wikipedia.org/wiki/FASTA_format
[2]: http://emboss.open-bio.org/
[3]: http://emboss.open-bio.org/html/use/
[4]: http://emboss.open-bio.org/html/use/apas01.html#SequenceFormatsFasta
Question: Thanks to the amazing package of KaryoploteR and a great tutorial website, I was able to make a rainfall plot.
https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html

However, one thing that I am suffering from making a perfect rainfall plot is color label is missing.
I would like to know and visualize what color is what like C>A is yellow, C>G is red for example.
Tutorial website states to ask a question to Biostars, I am posting my question here.

Thank you very much in advance.

Answer: You can use the standard `legend()` command as outlined in this issue here: https://support.bioconductor.org/p/124328/

Minimal example based on https://bernatgel.github.io/karyoploter_tutorial//Examples/Rainfall/Rainfall.html :

```
library(karyoploteR)

somatic.mutations <- read.table(file="ftp://ftp.sanger.ac.uk/pub/cancer/AlexandrovEtAl/somatic_mutation_data/Pancreas/Pancreas_raw_mutations_data.txt", header=FALSE, sep="\t", stringsAsFactors=FALSE)
somatic.mutations <- setNames(somatic.mutations, c("sample", "mut.type", "chr", "start", "end", "ref", "alt", "origin"))

somatic.mutations <- split(somatic.mutations, somatic.mutations$sample)
sm <- somatic.mutations[["APGI_1992"]]
sm.gr <- toGRanges(sm[,c("chr", "start", "end", "mut.type", "ref", "alt")])
seqlevelsStyle(sm.gr) <- "UCSC"

variant.colors <- getVariantsColors(sm.gr$ref, sm.gr$alt)
kp <- plotKaryotype(plot.type=4)
kpPlotRainfall(kp, data = sm.gr, col=variant.colors)
uniq_cols <- variant.colors[!duplicated(variant.colors)]
legend(x = "bottomright", fill=uniq_cols, legend=names(uniq_cols))
```

Now you can customise the call to legend() in any way you like.

Looks like this:![example rainbow plot][1]

[1]: /media/images/ab63d7b6-5198-4f0f-a350-d20d4ba4