cannin/biostars_qa huggingface dataset:
Question: I'am trying to open index.html from JBrowse after running ./setup.sh.
I have just loaded sample_data of volvox.
But some of the tracks are not loading file.
Showing "Error : HTTP 200 when fetching file"

Please help me in this
Thank You.!

Screenshot : https://photos.app.goo.gl/NQsa5SJC5w7APfHk8

Answer: First run
perl bin/flatfile-to-json.pl programme with features
and then run
run generatename.pl. program
Only after this results are displayed.
Question: I am attempting to run TransDecoder on transcriptome data generated by trinity. I have installed CD-HIT and exported the path correctly but TransDecoder fails when trying to run cd-hit-est.

Error, cmd: /data/xxx/src/cd-hit/cd-hit-v4.5.4-2011-03-07/cd-hit-est -r 1 -i /spin1

/users/xxx/trinity-transdecoder/redundant_top -o /spin1/users/xxx/trinity-transdecoder

/redundant_top.nr90 -M 0 -T 20 >/dev/null 2>/dev/null died with ret 256 at /usr/local

/apps/trinity/trinityrnaseq_r20131110/trinity-plugins/transdecoder/TransDecoder line

467.

However if I copy and paste this command and try and run it, it seems that cd-hit-est runs correctly and exits without error. I'm not sure what is wrong, if TransDecoder was unable to find the cd-hit-ext executable, it would crash before getting to this point, but it does not.

Answer: A perl exit code of 256 is the value from the wait call and should be shifted right by eight (C apparently does the same). So, an exit code of 256 corresponds to an exit code of 1, which means a general error. This could mean anything from errors in the code to permissions issues.

Try capturing the exit code in bash right after you run CD-HIT; it's possible it exits with a code you don't catch when running in the shell.

cd-hit-est .....
echo $?
Question: Hi all,
here are some questions about paired-end sequencing for NGS:

What are the main differences between mate-paired sequencing and paired-end sequencing; Should I care when I use tools like 'samtools', maq, etc.... ? Should one, and only one short read, should be paired with another one (1-1)?

What is removing duplicates ? does it mean that a pair of short reads has been mapped at two distint positions on the genome or does it mean that a pair matched too many time at one position ?

Knowing that bwa sampe "Generates alignments in the SAM format given paired-end reads. Repetitive read pairs will be placed randomly", is there any need to "remove the duplicates" ?

How does picard MarkDuplicates work ? How can I find the reads that have been 'tagged' ? will it remove the reads from the BAM file ?

Thanks

Pierre

Answer: Hm. Old questions, so nobody will read this, but I'm not entirely happy. Here are my answers:

Paired ends is supported by some technologies (Illumina and Sanger), where it is possible to sequence from both ends of a clone. Mate pairs involves making circular fragments using a linker sequence, and fragmenting them around the linker, and then sequencing the result. Illumina will read from each end of the fragment, 454 (and I believe Solid?) will read through it all.

Now this terminology isn't fixed, and lots of people will talk about mate pairs when doing paired end, or vice versa. Caution is advised! Also, I think Solid has some funky variations on this, but I haven't looked too closely. And, the mate pair protocol is rather unreliable, IME. Expect lots of non-mate-pair reads, and a wide range of insert sizes.

Oh, and I believe you can get mate pairs from 1.5K to 20K. We just got Illumina PE reads at 500bp inserts, this was considered experimental by the company doing it.

Removing duplicates can refer to several different things, but with all the second gen technologies, it is common to get a proportion of duplicate clones. This can skew things, e.g. for de novo assembly, it will give an artificially high coverage of a region, and it might be incorrectly identified as a repeat.

The random placement of sampe is probably to make sure that you get the right coverage for repetitive regions. So it's orthogonal to removing duplicates that are sequencing artifacts. And if you have a dinucleotide region of (AT), a read of (AT) and a read of (TA)* would not be considered duplicates, but be placed randomly.

Can't help you there, but Pierre seems to have it covered. :-)
Question: I have the tissue-specific RNA-seq data and its corresponding ATAC-seq data. I first identify all the differentially expressed genes. Then using ATAC-seq data to check whether the promoter regions of those genes are differentially accessible. It turns out that only a small portion of differentially expressed genes' promoter regions are differentially accessible. I try to find explanations. One reason I can think about is that the data is not perfect so missed some signal. What could be other reasons in the biology perspective?

Answer: We've never seen much in the way of correlation between chromatin openness at promoters and gene-expression changes. This is not totally unexpected. Studies of chromatin state find that chromatin in a promoter state is often the same across cell types (while enhancers tend to be more celltype specific).

It becomes obvious that gene expression is not entirely controlled at the chromatin openness level when one thinks about plasmid based reporter assays of promoter function. If you put e.g. AP-1 transcription factoer dependent promoter in a plamsid reporter, and then compare AP-1 activated vs non-activated, you'll see a difference in expression, despite the fact that plasmids don't have chromatin.
Question: Hi all,
I'm relatively new to bash scripting and working on one to map samples with BWA bash-loop.

My loop is the following:

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

./bwa mem -t 60 -c 1 -R '@RG\tID:${i}\tSM:${i}\tPL:illumina\tPU:Lane1\tLB:exome' MyReference.fq ${i}.fastq.gz ${i}.fastq.gz > ${i}.sam
done

The '-R' flag will add the read group info into the sam file. I am having trouble with having the actual sample name being added, rather than it just printed with a '${i}'. I've spending lot's of time googling, and trying to figure this out, but haven't been able to yet. Any advice would be appreciated.

Answer: Try `"` instead of `'` for quoting. Variable expansion needs double quotes.
Question: Hi Everyone,

I am trying to replicate the results of the following paper:
https://genome.cshlp.org/content/28/3/396

The data is long reads from Pacbio RS II and is available in SRA in the following location :
https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SAMN06328616&o=acc_s%3Aa

I am getting an error "ERROR, could not read quality metrics for FASTQ Sequence" when I try to use the program bax2bam to convert the legacy pacbio files to pacbio bam files.

Any suggestions how to convert the files from SRA to pacbio bam files to be used with the IsoSeq3 pipeline?

Thanks

Answer: As [already said][1] by [foeckingmf][2] it seems to be a bug.
I had the same problem and solved it downgrading bax2bam from version 0.0.11 to 0.0.9

conda install -c bioconda bax2bam=0.0.9

[1]: https://www.biostars.org/p/492792/#9463249
[2]: https://www.biostars.org/u/c0631eff/
Question: Good morning.

I am trying to use kmergenie v1.6972. But it is broken when I typed "make"

```
sudo make
g++ -o minia/Pool.o -c minia/Pool.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
make: g++: Command not found
make: *** [minia/Pool.o] Error 127
```

Firstly I install g++ like the next user told me. But it appears new errors. I have r-base, Rscript installed.

I show the error with `make`, `make check` and `./kmergenie`

```
make
g++ -o minia/Pool.o -c minia/Pool.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
g++ -o minia/Bank.o -c minia/Bank.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
In file included from minia/Bank.cpp:21:0:
minia/Bank.h:13:63: fatal error: zlib.h: No such file or directory
#include // Added by Pierre Peterlongo on 02/08/2012.
^
compilation terminated.
make: *** [minia/Bank.o] Error 1
```

```
./kmergenie /home/zhu/File.fq
running histogram estimation
Execution of specialk failed: [Errno 2] No such file or directory
fitting model to histograms to estimate best k
could not predict a best k value
Execution of decide failed (return code 0). If this is a fresh Kmergenie install, try running 'make check'.
```

```
make check
g++ -o minia/Bank.o -c minia/Bank.cpp -O4 -pthread -D_largeint -DKMER_PRECISION=4
In file included from minia/Bank.cpp:21:0:
minia/Bank.h:13:63: fatal error: zlib.h: No such file or directory
#include // Added by Pierre Peterlongo on 02/08/2012.
^
compilation terminated.
make: *** [minia/Bank.o] Error 1
```

Regards

Answer: sudo apt-get install g++
Question: How to do chromosome painting using [Satsuma synteny](http://satsuma.sourceforge.net/)?

Also how to download and visualize the result using mizbee?

Answer: Hope this may help you

SatsumaSynteny -t target_genome -q query_genome -n 2 –m 8 –o output_directory

I am also waiting for this command to finish and than I will explore Mizbee.
Question: I'm currently designing a workflow to analyze some exome sequencing data (VCF).

Now, I'd like to add OMIM data and HGMD data to my VCF file.

but I couldn't download some database that include them or find any software to annotate this data online or offline.

Do you know where I could find this kind of data file or which software doing those?

Regards.

Bahar

Answer: Hi Bahar,

We use the recommended tool SnpEff for variant annotation in VCFs, mainly to predict AA changes and effects. We also annotate the ID field with DbSNP using SnpSift (under the same project). However, to my knowledge there's also an option to annotate the VCF with several other clinical databases such as OMIM, check this link: http://snpeff.sourceforge.net/protocol.html

SnpEff: http://snpeff.sourceforge.net/

Hope this helps,
Alon
Question: Hi all,

I'm dealing with a fasta file with spaces at the end of line, which caused the problem. I didn't find a suitable way to remove them. Please kindly tell me the appropriate command for removing them?

Answer: Not a bioinformatics questions, you should try Stack Overflow for this, but here is a quick answer in perl:

perl -i.bak -pe 's/\h+$//' sequences.fa
Question: Hello,

I want to align multiple FASTA files in a directory with MAFFT.

Is it correct to do this with `parallel`(where `--thread n` is the maximum number of cores in the machine):

`ls *.fasta | parallel 'mafft --adjustdirection --thread n {} > {.}_mafft.fasta'`

If I am interpreting the line above correctly, each alignment will be performed in parallel using the maximum number of cores possible. Am I correct in setting the `--thread` option to the maximum number of cores for each alignment, or is `parallel` already taking care of that?

Or is it preferable to perform something like the for loop below and align the FASTA files sequentially:

for i in *.fasta; do
mafft --adjustdirection --thread n ${i} > ${i%.*}_mafft.fasta;
done

Thanks.

Answer: It would be more efficient to do the alignments sequentially using the maximum number of cores you have available (since you are able to use multiple threads). If you started multiple parallel jobs (each with multiple cores) they would compete for the same cores and the jobs would bog each other down.
Question: Some of the SNPs in my dataset are without rs ids, so in the bim file, they have dots instead of rs ids for them. So for errorless analysis how should replace the rs id column with that of the chromosome: location column instead?

Answer: This can also be done with plink 2.0's [--set-all-var-ids flag](https://www.cog-genomics.org/plink/2.0/data#set_all_var_ids).

But it's more valuable to understand the **awk** solution, since that opens the door to a broadly useful set of text-manipulation techniques.
Question: Hi,
I am analyzing T Cell Receptor sequencing data received from Adaptive Biotechnology which consists of summarized clonotype data which includes NT and AA sequence among many other variables. Some of the clones don't have an AA sequence but have a NT sequence. I believe these are called "unproductive" clones. I was wondering if anyone knows why some clones don't have an AA sequence (are unproductive) and more over how does Adaptive determine which clones don't have an AA sequence (what is their algorithm for this). Some of these clones that don't have an AA sequence have quite high frequencies, so I am particularly surprised that these don't have an AA sequence. One example of such as NT sequence is:

CAGAACTGGAGGATTCTGGAGTTTATTTCTGTGCCAGCAGCCTTTACAATGAGCAGTTCTTCGGGCCA

Thanks,
- Pankaj

Answer: This is a productive sequence, just load it to IgBlast and see for yourself:

![enter image description here][1]

' TRBV14 is productive and in frame with J

I think the situation can probably be explained as follows: the mapping software checks for conserved ``FGXG`` motif of Joining segment to ensure that J segment is mapped and as the last ``G`` is not in the read so the clonotype is not translated.

I've implemented a parser for ImmunoSEQ data a while ago in [VDJtools](https://github.com/mikessh/vdjtools): [executable](https://github.com/mikessh/vdjtools/releases/tag/1.1.2-SNAPSHOT), [docs](http://vdjtools-doc.readthedocs.io/en/latest/input.html#immunoseq). It converts it into a more compact format and checks if the rearrangements are truly non-productive, translating CDR3 sequences. I haven't tested it much so if anything goes wrong feel free to write in the issues section of the repository. Note that the option for conversion with ``VDJtools/Convert`` is either ``-S immunoseq`` or ``-S immunoseqv3`` depending on the version of Immunoseq software.

---

PS. It seems to be that it is a public clonotype: https://www.google.com/search?q=CASSLYNEQFF :)

PPS. Marked as Rasmussen encephalitis-specific clonotype http://www.nature.com/articles/ncomms11153 found in 0% control samples, possibly due to HLA restriction.

[1]: http://i.imgur.com/sib46Y8.png
Question: Hi all,

I have a weird question! I have some contiges and want to convert them to all possible illumina reads (suppose with length 100). I am wondering if there is any tool to do this?

Actually my problem is vise versa of the assembly problem. In assembly problem we have the reads and want to find the contigs. But here I have the contiges and want to find all possible small reads that can produce this contig.

Thank you,

Answer: While there are many tools out there, you can use `randomreads.sh` from [BBMap suite][1] with your contig. See the in-line program help for usage information.

[1]: https://sourceforge.net/projects/bbmap/
Question: Hi all,

I am interested in pulling out "CDS" fasta sequences from a GFF3 file. My GFF3 file lists one or more fasta sequences at the end.

I looked into using the `gffutils` python package but according to [this documentation][1], it does not deal with fasta sequences listed at the end of the GFF3 file.

Am I stuck writing my own parser, or is there a better solution?

Thanks for any suggestions!

EDITED to add a snippet of my GFF file:

```
##gff-version 3
##sequence-region Consensus_10_consensus_sequence_1 1 645439
Consensus_10_consensus_sequence_1	feature	gene	551	2104	.	+	.	name=gene
Consensus_10_consensus_sequence_1	feature	CDS	551	2104	.	+	0	name=YALI0C06512p CDS
...
Consensus_9_consensus_sequence_4	feature	mRNA	9625	9891	.	+	.	name=mRNA
Consensus_9_consensus_sequence_4	feature	CDS	9625	9891	.	+	0	name=YALI0A21351p CDS
##FASTA
>Consensus_10_consensus_sequence_1
TGCCACCTCCAAATTAACTCTCGCTTATTTCTTGTACCTGTCATATCACGTGATGTAGCT
TCCCAATCAAGAGCGGATCCTGCCTGTTTGGCTGCGTGGGTTTGCGTCTTCTTTCCGTTT
GAAGCAGTGGTATTATTCCCCCATTGTGCCAAAAGTAATGCTGAAAAGATGCCCACGAAT
>Consensus_10_consensus_sequence_2
CCGAAACCACAGCCATGATCAGAGTCACTCCTATTCCAACCCCGCCAACTGCCGTGGGGT
ACACGTTATATTGGGTAGGTGTGTATCCTTGAGACTTGAGCCAGGAGATCATGGAAGGTT
GGGTGCTGGCAATACAGGTGTTATTGTAGCAAAGAAAGATGAGAGAGAAGAAATAAATGT
GCCAGGTTTTCAAAGAGCGTTTCAAAACTTGCAAGAACGGCTCTTCTTTGTTACCGGTAT
CGCCGATTTGTGCTCTTCTCTCCTTGGCTATAGCAATGTCTTCCTTGGTGAAGTACCAGC
TGGTAGTTGTCTCCGGTGTGTTTGGGTTGACGAACATGGTGTAAAGTGCCACTGGAAATG
>Consensus_10_consensus_sequence_2
CACTCTCAAGCATTTAGGAACTTGTCAAGAGGTTCAAAGGTTGGAACTTGCAACTGAACT
GATCGCAACATAATCACCGCTATTAAACCCTGATTACAAGTGCTTCTGATTGCATCCACA
GTTCATTTCCATGGGCTAGGCTATACGAAAATACAAGGATTAGAAACTATATACAATTGA
CTCTGCAATCTTTCCCGCTAAACGGTGGTGTGGTTATGACCTGGCTCGTGTTCATGGCCG
...
```

[1]: http://pythonhosted.org/gffutils/examples.html#hybrid1-gff3

Answer: You should just split your file into two separate files, a GFF and FASTA. Then you can do:

https://gist.github.com/mdshw5/a8ce8748aa6f9e7a96607e7d8b2c8bfe
Question: I'd like to ask about the meaning of the 5th column (an integer number) in BED format results from repeatmasker.

9	100663131	100663387	LTR22_SS	475	+
9	70254161	70254685	ALTR2B_SSc	3460	+
9	96468811	96469391	LTR8_SSc	3756	-
9	116391614	116392469	LTR78	1152	+
9	4980341	4980930	LTR39_SSc	3300	-
9	16908116	16908359	MamGypLTR3	512	+
9	17432426	17432886	ALTR2B2_SSc	1914	-
9	18742941	18743430	LTR5_SS	2771	+
9	27131556	27132076	ERV3-1_SSc-LTR	969	-
9	30539515	30539909	LTR39B3_SSc	787	-

Searched around but did not find a definite answer, so seeking help here from repeatmasker experts. It is clearly not length of the repeat. Is it something indicating how repetitive this feature is in the genome? Or some identity scores indicating its match to some sort of repeat consensus?

Answer: It is the [Smith-Waterman alignment score][1] of matches with coding sequence.

I think (but am not 100% certain) a higher score means greater similarity with coding sequence, indicating a less repetitive region. It is used for cutoff filters, which are specific to different classes of repeats.

I don't know how useful it is to use these scores directly. Also see "[How to read the results][2]" from the [Repeatmasker documentation][3]:

> Smith-Waterman score of the match, usually complexity adjusted
The SW scores are not always directly comparable. Sometimes
the complexity adjustment has been turned off, and a variety of
scoring-matrices are used.

It may be worthwhile to [contact the developers directly][4].

[1]: https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/rmsk2bed.html#column-mapping
[2]: http://www.repeatmasker.org/webrepeatmaskerhelp.html#reading
[3]: http://www.repeatmasker.org/
[4]: http://www.repeatmasker.org/cgi-bin/form2mail?template=feedback.tmpl&title=Feedback%20Form
Question: I wonder if an annotation database exists which contains all known feature across the whole genom hg19 using keywords from [sequence ontology][1] ?

[1]: http://www.sequenceontology.org/

Answer: If you are interested in variants annotated according to the SO terms, download the [VEP file][1] from the [Ensembl FTP for GRCh37][2] (=hg19).

[1]: ftp://ftp.ensembl.org/pub/grch37/release-84/variation/VEP/
[2]: ftp://ftp.ensembl.org/pub/grch37/
Question: Hi,

I have one question, I need to get reads after filtration. Than normally I get health reads (health = reads after filtration, satisfy the conditions), but I need to get reads, which I remove.

For example I set filtration to get reads with mapQ bigger than 10, I used `samtools view -bg 10 ...` and I get reads with mapQ bigger than 10, but I want reads with mapq lower than ten.

Could you please help if it is possible. And other option, compare the raw bam file and filtrate bam file and get the difference between them.

Could you help me pls or do you know some tool to do this?

Answer: You can use bamtools to perform all sorts of BAM filtering tasks. Try the software at gkno.me (click on the 'Try it' button for simple instructions on getting the software) and you can get a while resource for next-generation sequencing analysis. To get the reads with mapQ less than 10, the following command line will work:

gkno bamtools-filter --in -mq "<10"

You can change the conditional command to easily get anything that you need. There are lots of additional filtering tasks (e.g. pulling out paired end reads, duplicate reads etc) that can easily be performed with gkno and bamtools.

As to your second question, what comparisons do you want to do? One easy thing to do, is apply the inverse of the above command:

gkno bamtools-filter --in -mq ">=10"

and you'll have all the reads not present in your filtered bam. You can get summary bam statistics using:

gkno bamtools-stats --in

and compare the results of the raw, and the two subsetted bam files.
Question: I want to use Salmon just to check the strandness of my samples. Is there a way to get the guess of libtype (--libType A) without getting all the output from the quantification? Or skipping running the quantification all-together?

Answer: https://salmon.readthedocs.io/en/latest/salmon.html#skipquant
Question: Hi,

There were already some discussion started to predict sample ethnicity from population genetic studies (1000G, gnomAD) e.g. https://www.biostars.org/p/300930/#300933 who used the PC from 1000G to create a predicitve model.

However I was wondering if there are already tools/R packages that take as input a vcf file and outputs the ethnicity of the samples ?

Thanks

Answer: Not R, but I would suggest taking a look at [peddy][1].

[1]: https://github.com/brentp/peddy
Question: I have 2300 vcf files from 2300 samples generated from haplotypeCaller of GATK from bacteria. Now i want to combine all VCF into one as multivcf sample file. I used CombineVariants from GATK but it is giving error when i give -V *.vcf.gz. I cant write -V file name 2300 times can anybody help me with this?

java -Xmx8G -jar /home/sark/GenomeAnalysisTK-3.8.1.0/GenomeAnalysisTK.jar -T CombineVariants -R ../ref/M._tuberculosis_H37Rv_2015-11-13.fasta --interval_padding 50 -V *raw.vcf.gz -L ../effluxgenes.bed -o multiple.vcf

Error:

##### ERROR
##### ERROR MESSAGE: Invalid argument value 'ERR040140_raw.vcf.gz' at position 8.
##### ERROR Invalid argument value 'ERR046796_raw.vcf.gz' at position 9.
##### ERROR Invalid argument value 'ERR046903_raw.vcf.gz' at position 10.
##### ERROR Invalid argument value 'ERR067581_raw.vcf.gz' at position 11.
##### ERROR Invalid argument value 'ERR067593_raw.vcf.gz' at position 12.
##### ERROR Invalid argument value 'ERR067606_raw.vcf.gz' at position 13.
##### ERROR Invalid argument value 'ERR067607_raw.vcf.gz' at position 14.
##### ERROR Invalid argument value 'ERR067608_raw.vcf.gz' at position 15.
##### ERROR Invalid argument value 'ERR067609_raw.vcf.gz' at position 16.
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
##### ERROR Invalid argument value 'SRR671854_raw.vcf.gz' at position 173.
##### ERROR Invalid argument value 'SRR671855_raw.vcf.gz' at position 174.
##### ERROR Invalid argument value 'SRR671866_raw.vcf.gz' at position 175.
##### ERROR ------------------------------------------------------------------------------------------

Can anyone help me with loop thing?

Answer: (...) -V *raw.vcf.gz (...)

GATK doesnt work like this.

You need

-V file1.aw.vcf.gz -V file2.aw.vcf.gz -V file3.aw.vcf.gz

**or** put all the paths in a file with the '.list' suffix

$ ls *raw.vcf.gz > vcf.list

$ java -Xmx8G -jar (...) -V vcf.list (...)