cannin/biostars_qa huggingface dataset:
Question: Dear all,
STAR manual said set --sjdbOverhang at readlength-1. I am confused by this with two questions.

1. If I generated genome index by setting it at 100 previously, and now I have a PE150 reads, should I re-generate genome index by setting it at --sjdbOverhang 149?

2. If I use --clip3pNbases 50 during mapping, should I set --sjdbOverhang at 99 (is from: 150 -50 -1)?

Thanks for your help.
Best.

Answer: You know the manual says about sjdbOverhang

> **In most cases, a generic value of 100 will work as well as the ideal**
> **value.**

http://labshare.cshl.edu/shares/gingeraslab/www-data/dobin/STAR/STAR.posix/doc/STARmanual.pdf
Question: Hello all,

I would like to remove any singletons/doubletons from my vcf file.

I can find all the singletons/doubletons using: `vcftools --vcf file.vcf --singletons --out file`

This output `file.singletons` which has this information:
`CHROM POS SINGLETON/DOUBLETON ALLELE INDV`

I thought I could extract the `POS` column from the `file.singletons`, remove the header and use as follow:

`vcftools --vcf file.vcf --exclude file.singletons --recode --recode-INFO-all --out final_vcf`

However this does not seem to remove any SNPs from the `.vcf` file. Is there anyway I can remove the SNPs listed in the `.singletons` file?

Thank you all for your help in advance

Answer: If anyone is having a similar problem - I found the solution!
Instead of using `--exclude` flag use the `--exclude-positions` flag for `vcftools`. Make sure you take the 1st and second column of the `--singletons` output.
Question: Hi all ~~

I'm qualifying my reads for further RNA-seq analysis.

When I did fastqc for the first time, I got failure in

- Per base sequence content
- Per sequence GC content
- Sequence duplication levels
- Overrepresented sequences
- Kmer content

Then I checked my over-represented sequences in blast and those over-represented sequences were for chloroplast genome. So I eliminated all chloroplast genome from my data.

Then I did fastqc again and I got failure in all I mentioned above except Over-represented sequences. Still I had 4 failures.
After I read through websites I understood failure in Sequence duplication levels is Ok because they might be due to highly expressed transcripts. So I felt happy for Sequence duplication levels.

On the other hand, the failure in Per base sequence content was due to first 13 bases. So I trimmed them and this module also healed.

But still I have failure in Per sequence GC content and Kmer content. Kmer content is even worse than before and the peaks are all over the positions while at the beginning peaks were around first 12 positions (which I trimmed them).

To make myself sure of lack of adapter sequences (The adapter content in Fastqc is completely alright) I put all adapter sequences used by Illumina in a file and tried to trimmed them off but as I knew there were no adapters and nothing changed in quality.

So do you have any suggestion for having more qualified reads? Did I do something wrong? Why Kmer became worse?

Answer: For RNAseq datasets one typically sees failures in those modules. I should point out that FastQC is really geared toward whole-genome sequencing. We use it for all types of datasets, but many of them will have failures in various modules that can be ignored. In particular, the modules you mentioned will commonly fail on RNAseq samples and I wouldn't worry about that. Additionally, I would encourage you to not trim off the first 13 bases that are biased. The bases are actually correct and you'll lower your mapping rate (and cause some false alignments) by doing that.

BTW, remember with RNAseq data you expect anything short and highly expressed to be sequenced many many many times. That alone will cause a shift in the GC and kmer profile and cause "failures" in these tests. Don't worry about stuff like that.
Question: I have a .gtf file as below I would like to add "chr" to the first column of the file but not in first 5 rows?

#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
1 ensembl chromosome 1 300239041 . . . ID=1;Name=chromosome:AGPv1:1:1:300239041:1
1 ensembl exon 3 104 . + . Parent=GRMZM2G060082_T01;Name=GRMZM2G060082_E07

output:

#!genome-build GRCh37.p13
#!genome-version GRCh37
#!genome-date 2009-02
#!genome-build-accession NCBI:GCA_000001405.14
#!genebuild-last-updated 2013-09
chr1 ensembl chromosome 1 300239041 . . . ID=1;Name=chromosome:AGPv1:1:1:300239041:1
chr1 ensembl exon 3 104 . + . Parent=GRMZM2G060082_T01;Name=GRMZM2G060082_E07

I used the foloww cods but it add "chr" to the first 5 lines, as well.

cat Homo_sapiens.GRCh37.gtf | sed 's/^/chr/' > chr.gtf

Answer: awk 'OFS="\t" {if (NR > 5) $1="chr"$1; print}' in.gtf

Could have been found on google easily...
Question: I'm using UpSetR to create UpSet diagrams [like this][1], but I can't find a way to add a main title to the plot. Does anyone know if this functionality is possible in the UpSetR package at the moment?

[1]: http://caleydo.org/assets/images/projects/upset//upsetr-banana.png

Answer: Half-baked solution is to add the title to the image itself using the [convert][1] tool part of the ImageMagick suite. At least you can do it programmatically. For example:

convert -gravity North -annotate +0+0 'My Title' -pointsize 20 input.png output.png

This will add the text "My Title" on the upper side of the image (`-gravity North`) right in the middle. It can handle several image formats in input and output. ImageMagick is useful suite to have anyway...

**EDIT**

To choose a font just add the option `-font`, for example `-font TimesNewRoman`. To see available fonts you can use `convert -list type`.

[1]: https://www.imagemagick.org/script/convert.php
Question: Hi.

I made vcf files using UnifiedGenotyper in GATK and found that many of mutants had 1/2 at genotype (GT) field.

As I understand, 0/0 = homo ref, 0/1 = hetero, 1/1 = homo alt.

Does anyone know what exactly 1/2 mean?

Thank you!

Answer: just to expand Devon's answer, which is perfectly correct...

this is due to a VCF format constrain, which states that 0 must always be the reference allele. this doesn't mean that the reference allele must be in a particular sample, it just means that a particular base is used as reference by a scientific community agreement.

in case you're working with human samples, this agreement is now supported by the Genome Reference Consortium, and I really hope that when the time comes (we shouldn't be that far) they'll be able to build a consensus reference based on entire world populations allele frequencies, but meanwhile we have to stick to the current reference. current GRCh38 [has already done][1] some work in this direction though.

since humans are diploids, sequencing a particular position will give you a pair of alleles, and usually the reference will be 1 of that 2 alleles. but this is not always the case, and you may find (as you already did) that a particular position doesn't contain the reference allele. this is rare, but far from being impossible due to the lack of population genetics evidences in the reference selected, as explained previously.

[1]: https://www.biostars.org/p/94826/
Question: How to convert bam to fasta?

I have bam files and reference fasta. I want to get fasta sequence from bam file:
![see picture][1]

Is there a software that can do this?

[1]: http://s22.postimg.org/hwdzkaktt/fasta.png

Answer: What the original poster needs is creating a consensus sequence out of BAM alignments. It is easy to misread that into the problem of extracting the read sequences from the BAM file.

See this for an answer: https://www.biostars.org/p/184944/#185560
(notably the title of this post is incorrect as well, the OP there wants FASTA file not FASTQ)

But note that the problem can be more complicated than what it initially seems. It only works in ideal cases.

In the most general case the SNPs may not resolve into a single consensus sequence because that is not how the SNP callers work. They optimize the SNPs for a local region and not against a consensus.
In addition in many (most) cases you may have a population of cells hence SNPs may represent different DNA molecules.

This is not to say that there is no value in producing a consensus but it is good to keep in mind the limitations of the the process.
Question: Hello,

I have a problem using ExomeDepth with a referenceFASTA file. I analysed app. 500 target genes (>7000 Exons, hg19) and without including a FASTA file and the GC content everything worked fine.

As GC content influencing amplification, coverage and CNV detection, I want to include this factor. Therefore I tried several things to get appropiate FASTA sequences.

BioMart

mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl")
FASTA <- getSequence(chromosom,start,stop,type="hgnc_symbol",seqType="gene_exon", mart = mart)

Error:

`Reference fasta file provided so exomeDepth will compute the GC content in each window
Error in (function (classes, fdef, mtable) :
unable to find an inherited method for function ‘scanFa’ for signature ‘"data.frame", "GRanges"’`

BSgenome

FASTA<-getSeq(BSgenome.Hsapiens.UCSC.hg19,chromosom,start,stop)

show(FASTA)

A DNAStringSet instance of length 10 width seq
[1] 845 GTTGAAAAGTGATCAGGTTCATTTTATTGACTACACAGAAGCAATTCCATTT...GAGGAGGCAGATCACGGCGAAGACAATGAAGCTGTACGGGCCGAGGCCCTC
[2] 129 CCTGGATGAACGGGAAGATCAAGCCCACGGTGAAGTTGGAGAGCCAGTGCAC...GCCGAGAGGACTGCAGGAAGATCTCAGTGATGAGCAGCGCGGGTATGGGAC
[3] 77 CTGGGCCCGAGGGCATGTCCTATGACGTAGGAGATGACACAGACGATGCTGATGTATGGCATCCAGGACACTGTGTC
[4] 103 CCTGCAGTGCCAGAGCTGCAGTGAGCACGCAGCAGGCTATGAGGCAGATGGAGAAGCCCAGCAGCAGCAGCAGCCTCCGACCCAGGAGCTCCACCACGAACAC
[5] 112 CGGCGCAGAAGGTCATGACCACGTTCACGGCCCCGGTGCCGGCCGTCACGTA...CTCCTCCGGCACGCCGGCGCTCAGGTAGATCTGGTCCGCGTAGTAGTAGAT

Error:

Reference fasta file provided so exomeDepth will compute the GC content in each window
Error in (function (classes, fdef, mtable) : unable to find an inherited method for function ‘scanFa’ for signature ‘"DNAStringSet", "GRanges"’

3. Downloading human_g1k_v37.fasta.gz

Chromsome, Start & Stop including only the 7000 exon locations.
With method 1 and 2, I was able to get target sequences, but the counting (getbamcounts) didn't work, different errors occur, warning something is wrong with the FASTA file. I think there are some file formatting issues.

my.counts<- getBamCount(data.frame(Chromsome, Start,Stop),
bam.files = bam.files,
include.chr = F ,
referenceFasta =FASTA)

I dindn't tried the human_g1k_v37.fasta.gz so far, because I don't know how to load it in R.

Do you have an idea how to transform the files that they work or how to load the whole genome FASTA File?

Answer: The reference FASTA file should be the same as that used for your mapping step, so you should have it already.
Question: Hi, how do you map reads using TopHat when you have several replicates? Does it work like this:

```
tophat -p 8 -o output_dir \
./path_to_bowtieindex \
./path_to_replicate1 \
./path_to_replicate2 \
./path_to_replicate3
```

I am omitting most of parameters I use to keep it short.

Answer: Technical or biological replicates? What you showed will work for the former, but not the latter. For the latter, you need individual runs of tophat, since you need the output in different files.
Question: I have a very large text file (9.2GB) that contains data in two columns (named All_SNPs.txt). I would like to use this file to convert my chr:pos to rs in my plink files. The file looks like this:

1:116342	rs1000277323
1:173516	rs1000447106
1:168592	rs1000479828
1:102498	rs1000493007

However, `plink` produced an error stating that All_SNPs.txt contains duplicates in column 1. I have tried different ways of removing the duplicates; specifically, I attempted the following line commands:

awk '!x[$1]++ { print $1,$2 }' All_SNPs.txt > All_SNPs_nodup.txt

sort All_SNPs.txt | uniq -u > All_SNPs_nodup.txt

cat All_SNPs.txt | sort | uniq -u > All_SNPs_nodup.txt

However, each time I faced the same two problems: 1) either, the code will not make any difference and I will receive the same file back (this is for line commands that uses `cat`) or 2) I will receive an error : the procedures halted: have exceeded (this is for `sort | uniq`, and `awk`).

I will be very grateful for any ideas of how I can make this work. Thank you very much.

PS, this file is far too large to open it in R.

Answer: Split the text file into smaller one using the `split` command. You could split by size, as example:

split -b 200m filename

This will produce files named ' xaa, xab, xac..' . Now use awk, but with simpler syntax

awk -F"\t" '!seen[$1]++' xa*

And after that, join files using a sample `cat` into the destination file.
Question: **Where can I download the population frequency data of all mitochondrial variants?**

I checked [the HapMap database][1], but it only includes the nuclear variants.

[1]: http://hapmap.ncbi.nlm.nih.gov/downloads/frequencies/latest_phaseII_ncbi_b36/fwd_strand/non-redundant/

Answer: The GenBank frequency data (out of about 30k full length mitochondrial sequences) for any given variant is available in [MITOMASTER][1]

MITOMAP provides curated polymorphisms and disease dumps of various formats with these counts and frequencies

- http://mitomap.org/cgi-bin/disease.cgi?format=vcf
- http://mitomap.org/cgi-bin/polymorphisms.cgi?format=vcf

These will probably evolve as we develop breakdowns by major haplogroup.

[1]: http://mitomaster.mitomap.org
Question: Heys,

I'm working with two ssh clusters and I want to transfer some data from one to another. I could download it to my pc and then to the other cluster, but is more than 100 GB and I imagined there would be an easy way to transfer the data between both clusters.

I've been trying to look for info in internet, but I didn't fully got it. If anybody could help me would be great!

Thanks in advance!

Answer: Did you checked this:
https://unix.stackexchange.com/questions/106480/how-to-copy-files-from-one-machine-to-another-using-ssh
Question: I have been downloading a lot of NCBI short read archive files (.sra files) and I know how to convert them to bam, sam, fastq, etc but is there a way to generate my own .sra files? The tools in the SRA toolkit do not seem to convert to .sra

​Relevant Information:

- [SRA toolkit][1]
- [Handbook][2]

[1]: http://eutils.ncbi.nih.gov/Traces/sra/?view=software
[2]: http://www.ncbi.nlm.nih.gov/books/NBK47540/#SRA_Download_Guid_B.5_Converting_SRA_for

Answer: I have never used it but Mothur has a command "make.sra".
Question: I have a bash script about HaplotypeCaller.sh. Is it possible to run this in gnu parallel ?

cat HaplotypeCaller.sh | parallel

#!/bin/bash

for i in *.bam; do
gatk --java-options "-Xmx40g" HaplotypeCaller -R /media/gatk/Homo_sapiens_assembly38.fasta -I $i -O $i.vcf.gz --dbsnp /media/gatk/dbsnp_138.hg38.vcf.gz -L /media/gatk/target/ag_V6_list.interval_list; done

Answer: nextflow (not tested)

params.ref="/media/gatk/Homo_sapiens_assembly38.fasta"
params.bams=""
params.dbsnp="/media/gatk/dbsnp_138.hg38.vcf.gz"
params.intervals="/media/gatk/target/ag_V6_list.interval_list"

Channel.fromPath(params.bams).splitCsv(header: false,sep:',',strip:true).map{T->file(T[0])}.set{bamfiles}

process hapCaller {
tag "${bam.name}"
memory "40g"
input:
file bam from bamfiles
output:
file("${bam.getBaseName()}.vcf.gz") into vcf
script:
"""
gatk --java-options " -Xmx${task.memory.giga}g -Djava.io.tmpdir=." HaplotypeCaller \
-R ${params.ref} \
-I {bam.toRealPath()} \
-O "${bam.getBaseName()}.vcf.gz" \
--dbsnp "${params.dbsnp}" \
-L "${params.intervals}"
"""
}
Question: Hi,

Having two BAM files from NGS data, how can one check if they are the BAM files (left and right) from a paired end mapping of the same sample?
Thanks for the help.

Answer: Mates in a pair shouldn't be mapped separately. To check if a BAM file has reads that arise from paired-end sequencing (and were mapped accordingly), just check for the 0x1 bit in the FLAG field (i.e., see if samtools view -f 0x1 something.bam prints any reads).
Question: Hi,

I have a file (file1.txt):

AHR_si	liver
AHR_si	liver
AHR_si	liver
AHR_si	large_intestine
AHR_si	liver
AHR_si	large_intestine
AHR_si	liver
AHR_si	skin
AHR_si	liver
AHR_si	pancreas
then the file continues.......

I want to count the number of occurrences of cervix (not shown here as the file1.txt is showing as head -n10) appears against each column1. So:

grep cervix file1.txt | uniq -c

1 AIRE_f2	cervix
1 ARI3A_do	cervix
1 FOSB_f1	cervix
1 FOXQ1_f1	cervix
1 HEN1_si	cervix
1 HNF4G_f1	cervix
1 JUNB_f1	cervix
1 NFAC1_do	cervix
1 NR2F6_f1	cervix
1 PTF1A_f1	cervix
1 ZN350_f1	cervix

The above is the total output. As you can see there is not even a single occurrence of AHR_si with cervix. But I still I want the output like this:

0 AHR_si cervix
1 AIRE_f2	cervix
1 ARI3A_do	cervix
1 FOSB_f1	cervix
1 FOXQ1_f1	cervix
1 HEN1_si	cervix
1 HNF4G_f1	cervix
1 JUNB_f1	cervix
1 NFAC1_do	cervix
1 NR2F6_f1	cervix
1 PTF1A_f1	cervix
1 ZN350_f1	cervix

There are many other column1 value where cervix did not matched. So, I want in my output those lines as well with occurrence of zero.

Thanks in advance,

Waqas.

Answer: In bash you would have to use gawk, setting an array to count the occurrency of each sample name.

But for complicated tasks I would prefer to use R:

> install.packages("dplyr")
> library(dplyr)
> f = read.table("file1.txt")
> f
V1 V2
1 AHR_si liver
2 AHR_si liver
3 AHR_si liver
4 AHR_si large_intestine
5 AHR_si liver
6 AHR_si large_intestine
7 AHR_si liver
8 AHR_si skin
9 AHR_si liver
10 AHR_si pancreas
11 AIRE_f2 cervix
12 ARI3A_do cervix
13 FOSB_f1 cervix
14 FOXQ1_f1 cervix
15 HEN1_si cervix
16 HNF4G_f1 cervix
17 JUNB_f1 cervix
18 NFAC1_do cervix
19 NR2F6_f1 cervix
20 PTF1A_f1 cervix
21 ZN350_f1 cervix

> f %>% group_by(V1) %>% summarise(n.cervix=sum(grepl('cervix', V2)), n.liver=sum(grepl('liver', V2)))
# A tibble: 12 x 3
V1 n.cervix n.liver

1 AHR_si 0 6
2 AIRE_f2 1 0
3 ARI3A_do 1 0
4 FOSB_f1 1 0
5 FOXQ1_f1 1 0
6 HEN1_si 1 0
7 HNF4G_f1 1 0
8 JUNB_f1 1 0
9 NFAC1_do 1 0
10 NR2F6_f1 1 0
11 PTF1A_f1 1 0
12 ZN350_f1 1 0
Question: I think I found a bug in the R package of "ChAMP". Even In its latest version, this bug were still existed there. As the following code shown in line 237 of function `champ.lasso`,

myDf <- data.frame(myDf, dmr.beta.means[,22:24]) #

Actually, I think it should be

myDf <- data.frame(myDf, dmr.beta.means[,17:19]) #

That's my opinion, I have contact the authors, However, there is no any response for long time. Is there any other friend found this problem/bug?

A: Package can be found here: http://www.bioconductor.org/packages/devel/bioc/html/ChAMP.html

Answer: Hi Shicheng,

Thank you for your question. This update has been submitted and should be available from Bioconductor by the end of the week.

Best wishes
Tiffany
Question: I am trying to pipe a QuerySequences using bash to perform a blastn on a ReferenceGenome. The problem is that I have to many sequences and, to avoid producing a file with to unnecessary registers (memory issues), I want to limit my output to those results with a percent of identity of 90% or higher.

I have run the following script without the flag I am going to mention afterwards and it works perfectly.

#-perc_identity 90?????
blastall -p blastn -d ReferenceGenome -i QuerySequences -G 1 -E 2 -W 15 -F "m D" -U -e 1e-20 -m 8 -a 8 -o NAME.blast.out

But when I try with -perc_identity 90 flag the following error appears:

[blastall] ERROR: Arguments must start with '-' (the offending argument #4 was: '90')

I proved the flag in several position (after the blastn, after Querygenome, after "m D") and the only thing that changes is the number after the # in the error message.

Does anybody knows the probable reason of this error and how to solve it?

Thank you very much for your help.

Answer: You are mixing BLAST versions. "-perc_identity" works for the newer BLAST+, but you are using legacy BLAST. I don't think legacy BLAST has an option for percent identity, so if you want it, you should download and install BLAST+.

edit: be aware you'll have to change your whole command-line.
Question: Hi everybody,
I want to convert my bedgraph (from MACS2) into a bigwig, is there any available program different than bedgraphtobigwig (I have problems with the instalation)

Thanks

Answer: There should not be any installation. You just download the binary and run it.
Question: Hi,

I have problem with installing biomaRt for R version 3.0.2 in windows. does anyone has idea what would I do?

Installing package into ‘C:/Users/Sky/Documents/R/win-library/3.0’
(as ‘lib’ is unspecified)
Warning message:
package ‘biomaRt’ is not available (for R version 3.0.2)

Answer: You are using biocLite() to install and you have run the following command?

source("http://bioconductor.org/biocLite.R")
Question: Hi

I have two `transcripts.gtf` files which are outputs of cufflinks. One is 69.4 Mb and the other one is 38 Mb.

When I merge these files with their annotation file the `merged.gtf` file is 75 Mb. Why it is small in size? (69.4 + 38 = 107.4), so I mean the minimum size of merged file should be 107.4. Isn't it?

Answer: did you merge them with cufflinks? if so, the transcripts that are common to both files will be present only once in the merged.gtf. that's what merging is for. also, counting the number of transcripts is more informative than looking at the size.
Question: I have some large bam files for a pipeline. But the RAM ran out when I was performing the pipe. Sorting bam will cause fatal bugs to my pipeline and cannot be debuged for now, hence not the solution. The only option is to downsize the input (upgrading the RAM cannot be done for now due to the high expense).

Then I thought I could split the original bam to several smaller bam files. But I cannot skip the sorting step. E.g. `samtools view` requires sorting and indexing before spliting.

Is there other way to do this?

Thank you.

Answer: Using samtools *alone* the bam file needs to be sorted. However at the cost of materializing the actual string values, you can use awk to do this:

```
samtools view -h $bam | awk 'substr($1,0,1) == "@" || $3 == "chr1"' | samtools view -hb > chr1.bam
```
etc.
Question: Is there other pathway analysis tools except GSEA using gene set from MSigDB database? For example, tools can use user-defined pathway gene set or support gene set (gmt) from MSigDB directly? Thank you.

Answer: You could use `qusage` http://clip.med.yale.edu/qusage/. This is a Bioconductor package and is easy to use with R.

You could also use the `piano`: http://www.bioconductor.org/packages/release/bioc/html/piano.html. This is also a Bioconductor package and easy to use.

Go through the manuals and try to follow steps with the given examples.
Question: I just received my genotyping report files from Illumina, with the locus summary file, and I want to convert the files to Plink files.

They came as CSV so I converted them to tab delimited (.txt) and wanted to use gcta to convert into map + ped. However, I get the map file but get an error ("segmentation error") when it tries to create the .ped

**My locus summary file**

```
Index Name Chr Position AA Freq AB Freq BB Freq Call Freq Minor Freq 10% GC 50% GC
1 1:10002775-GA 1 10002775 0.00 0.00 1 1 0.00 0.4114687 0.4114687
2 1:100152282-CT 1 100152282 0.00 0.00 1 1 0.00 0.3868383 0.3868383
3 1:100154376-GA 1 100154376 0.00 0.00 1 1 0.00 0.4968892 0.4968892
4 1:100154844-CA 1 100154844 0.00 0.00 1 1 0.00 0.3769583 0.3769583
```

**Genotyping report file**

```
[Header]
GSGT Version 1.9.4
Processing Date 3/12/2015 2:17 PM
Content MEGA_Consortium_15063755_B1.bpm
Num SNPs 1548495
Total SNPs 1705969
Num Samples 572
Total Samples 576
File 1 of 6
[Data]
SNP Name Sample ID Allele1 - Forward Allele2 - Forward GC Score X Y
1:10002775-GA WG0238334-DNAA01_Control G G 0.4115 0.027 1.531
1:100152282-CT WG0238334-DNAA01_Control C C 0.3868 0.072 1.651
1:100154376-GA WG0238334-DNAA01_Control G G 0.4969 0.019 0.855
```

Is there any way of doing this? I would appreciate any suggestions.

Thank you

Answer: Try this

- https://www.biostars.org/p/13302/
- https://www.biostars.org/p/51928/