cannin/biostars_qa huggingface dataset:
Question: I would like to double check whether to use VCF concatenate or VCF merge on my chromosome files. I have done SNP calling using FreeBayes, but split this by chromosomes in order to call SNPs in parallel. I also split some particularly large chromosomes by chromosome position with no overlap, e.g. 1:500,000; 500,001:1,000,000, with the end result being:

`Chromosome1-position-1:500,000`

`Chromosome1-position-500,001:1,000,000`

`Chromosome2`

`Chromosome3`

I want to first combine the separate VCF files for single, extra data heavy chromosomes that have been split into multiple VCF files for processing. I.e.

`Chromosome1-position-1:500,000` + `Chromosome1-position-500,001:1,000,000` --> `Chromosome1`

Then I want to combine all of the separate VCF files into 1, i.e.

`Chromosome 1 {merged from step 1}` + `Chromosome 2` + `Chromosome 3` --> `SNPs`

The VCF tools manual leads me to believe that using `VCF-concatenate` is the appropriate command for both as they are all separate files of separate chromosomes that I just need to re-attach back to each other but I'm unsure if this is the case. Any advice would be appreciated.

Answer: Take a look at `bcftools concat` (https://samtools.github.io/bcftools/bcftools.html#concat).

In general, concatenate means adding on rows to your vcf (e.g. re-combining split chromosomes), while merge means creating a superset of variant calls across multiple individuals. Concatenating is relatively straightforward, you just need to keep your headers straight (e.g. going from multiple files each with their own headers, to one file with one set of headers, and ideally an additional header row documenting the command used to concatenate).

Merging across individuals can be more problematic. If you have a gvcf, which represents all positions whether they contain a variant or not, you can merge fairly easily. However, if you only have vcfs, which only report a genomic location if there is a variant in that individual, you encounter a missing data problem. When a variant is reported in A.vcf, but not in B.vcf, the merged file will record the variant as missing "./." for sample B. Does that mean there was insufficient coverage to make a call, or was there plenty of coverage and simply no variant reads? It doesn't sound like this is what you're doing here, so don't worry about it for now. :)
Question: I am doing some quality control on a bunch of affymetrix chip and I noticed a consecutively bright spot in every single one of these chips. I assume this has a technical use since I observed this in different chips. Can anyone tell me what it is? The image below is an example with exceptionally low overall intensity but that spot in the middle is shiny as always. Thus particular chip is Affymetrix U133 Plus 2.0. The image is generated by `image()` function of R on the affy object.

Answer: That's the control probe set. You **should **see that on every chip, so it's useful as a QC indicator.
Question: I am new to Hi-C analysis.I am working with a test data set.can any one give me a clear work flow of Hi-C analysis and also can any one tell me about the maths behind the expected matrix and inter and intra chromosomal maps

Answer: You're essentially asking "can someone tell me everything I need to do and why", while also not saying what you are aiming to do with this data.

Have you read the key papers? Yaffe and Tanay (2011), Lieberman Aiden et al. (2009), Dixon et al. (2012)? The supplementary of Dixon et al. is quite detailed. There's also a decent quick tutorial hosted at BioConductor (pdf, R) that walks you through extracting PC eigenvectors from Lieberman's intrachromosomal data and getting some matching chromatin features.
Question: So, I am trying to split a SAM file into many smaller SAM files -- by line number, where the line number is a multiple of 8.

I tried the following series of commands:

samtools view -H Mapped_Pila/LBG07/LBG07_pe.sam > Mapped_Pila/LBG07/LBG07_pe.header.sam && samtools view Mapped_Pila/LBG07/LBG07_pe.sam | split -l 8000000 - Mapped_Pila_Split/LBG07/LBG07_pe.split. && find -L . | grep Mapped_Pila_Split/LBG07/LBG07_pe.split. | parallel --gnu -j4 "echo Mapped_Pila/LBG07/LBG07_pe.header.sam {} >> {}.tmp && rm -f {}" && rename "s/\.tmp$/\.sam/" Mapped_Pila_Split/LBG07/LBG07_pe.split.*

Unpacking this:

1) get the header in a separate file

2) get the headerless SAM alignments and split them by a multiple of 8 (8000000)

3) find the split files from the original files, and for each of them, cat the header, then the headerless records into a new file

4) Remove the original split file

5) rename the tmp files to sam

This results in a SAM file that is parsed as truncated by samtools. So, I think I am missing something.

Generally, I want to split sam files that are very large into more tractable parts and then later recombine them with samtools merge.

Answer: You should add the header back to each of the splitted files.

something line `cat Mapped_Pila/LBG07/LBG07_pe.header.sam split_file > split_file.sam`

It should be

cat Mapped_Pila/LBG07/LBG07_pe.header.sam {} > {}.tmp && rm -f {}
in your command.
Question: I’m searching for a downloadable annotation file for different types of ncRNA: in particular, rRNAs and retrotransposons for mouse. Most ncRNAs are not a problem – for instance, Ensembl provides a FASTA download of ncRNAs. Luckily, the actual sequence data is rather short and the file contains full annotation.

Unfortunately, it doesn’t contain retrotransposons. I’m particularly interested in identifying LINEs and SINEs. Where could I get that information from? I tried Ensembl Biomart but was unable to find the appropriate filter there. The Ensembl API also contains a promising method – TranscriptAdaptor::fetch_all_by_biotype but I cannot find a list of biotypes anywhere, and if I read the ncRNA documentation correctly, Ensembl simply doesn’t annotate retrotransposons at all.

Is there any way to obtain that information?

Answer: Answering my own question:

Ensembl Biomart actually allows filtering by Gene type: “retrotransposed”.

However, this only works in the current build of Ensembl and I am unfortunately forced to use an archived build since I’m working on the genome version NCBIM37.
Question: Hello,

I'm trying to align some test .fastq files to a reference using kallisto.

The reference is composed of human transcriptome + a couple of plasmid sequences (~12000 characters) stored in .fa format. I generated the index using this command:

humFa="/path/to/ucsc/fasta/files/RefGenomes/H_sapiens/hg19/*.fa"
plasFa="/path/to/plasmid/fasta/files/*.fa"

kallisto index -i humPlas_kallisto_transcripts.idx $humFa $plasFa --make-unique

...and the resulting file is. 70.49 GB.

I have tried to align paired end .fastq files to this index using kallisto, but I keep running into issues:

# On my Mac laptop (macOS 10.13.3, 3.5 GHz Processor, 16 GB memory):

The issue seems to be that the program dies prematurely, but I don't know why. I run this with only Terminal open, and I don't touch anything while it's running:

./kallisto quant -i ~/Desktop/humPlas_kallisto_transcripts.idx -o ~/Desktop/kallOutput/ -b 100 ~/Desktop/tstFastq/R1.fastq ~/Desktop/tstFastqR2.fastq

[quant] fragment length distribution will be estimated from the data
[index] k-mer length: 31
[index] number of targets: 50,798
[index] number of k-mers: 2,969,625,638
Killed: 9

# On a computational cluster:

I run the same command (as above):

kallisto quant -i $kallistoIdx -o $outputFileLoc -t 4 -b 100 $Read1 $Read2

with the following cluster queue settings (SLURM):

#SBATCH --nodes=4
#SBATCH --ntasks=4
#SBATCH --cpus-per-task=4
#SBATCH --mem=256G

...and then the job auto-aborts after six hours because it hasn't completed in that time.

Am I correct in thinking this kallisto alignment is taking suspiciously long/is being suspiciously buggy? Has anyone run into either of these issues? Am I missing anything that might be making kallisto slower?

Thank you for your help!

# EDIT - SOLUTION
Thanks for the feedback - indeed, the problem is that I was using UCSC's *genome* files, not *transcriptome*.

I got the transcriptome corresponding to hg19 from the Ensembl archive here:
wget ftp://ftp.ensembl.org/pub/release-67/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh37.67.cdna.all.fa.gz

Then I regenerated the index using that file in my $humFa path. The resulting index was much smaller.

Now everything is working well, and I can align to my plasmids+human genome. Appreciate the help!

Answer: > The reference is composed of human **genome** + a couple of plasmid sequences (~12000 characters) stored in .fa format.

kallisto uses a reference **transcriptome**.
Question: Hello everyone,

I have a question regarding my single-cell RNA-seq data. I have the following pair-end data in `fastq.gz` format.

Read1 (contains 6bp UMI, followed by 6bp cell barcode info and the rest is a polyT stretch):

@J00182:79:HV2WWBBXX:6:1101:11160:38873 1:N:0:ACAGTG
GAGAAGACAGTGTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTCTT
+
AAAFFJJJJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJ--A-----7AJJFAF-AJJJJJJJJJJJ

Answer: Obviously as the UMI-tools author, I'm going to agree with @Devon Ryan! However, if you really want to do it this way I can think of two options:

**Method 1: Extract the the UMI and CB using UMI tools and then reformat the headers. So, using the Biopython SeqIO fastq parser (I'm sure better parsers are available)**

from Biopython import SeqIO
outfile = open("output.fastq", "w")
for read in SeqIO.parse("my_extracted_read.fastq", "fastq"):
header = read.id
read_name = header.split()[0]
id, cb, umi = read_name.split("_")
new_name = "%(cb)s_%(umi)s#%(id)s" % locals()
read.id = new_name + " " + " ".join(read_name[1:])
SeqIO.write(read, outfile, "fastq")

outfile.close()

There are definaly more efficient ways to write this, but I think this is the clearest. I also don't remember how `SeqIO` deals with the `@`, so you'll need to look into that.

**Method 2: Write your own parsing code.**

Take a look at `extract.py` from `umi_tools`, its not particulalry complex, and it should give you an idea for how to do the processing yourself.
Question: I have a large vcf.gz file for 70 samples. It contains many positions where all 70 samples match the reference. I'd like to filter those lines out.

I tried:

bcftools filter -e '%ALT == .' ./SIDXXX.vcf.gz

but it came back with ...

[filter.c:1741 filters_init1] Error: the tag "INFO/%ALT" is not defined in the VCF header

Thanks!

Answer: try to use the variable `N_ALT` see https://samtools.github.io/bcftools/bcftools.html (with N_ALT==0)

> variables calculated on the fly if not present: number of alternate
> alleles; number of samples; count of alternate alleles; minor allele
> count (similar to AC but is always smaller than 0.5); frequency of
> alternate alleles (AF=AC/AN); frequency of minor alleles (MAF=MAC/AN);
> number of alleles in called genotypes; number of samples with missing
> genotype; fraction of samples with missing genotype
>
> N_ALT, N_SAMPLES, AC, MAC, AF, MAF, AN, N_MISSING, F_MISSING
Question: Hello
i have a bed file like this

1	3004104	3004110
1	3005819	3005825
1	3008315	3008321
1	3008893	3008899
1	3009812	3009818
1	3012422	3012428
1	3015794	3015800
1	3016183	3016189
1	3024019	3024025
1	3025279	3025285

and i would like to create another one out of this to in this format :

1	1	3004110
1	3004104	3005825
1	3005819	3008321
1	3008315	3008899
1	3008893	3009818
1	3009812	3012428
1	3012422	3015800
1	3015794	3016189
1	3016183	3024025
1	3024019	3025285

Do you know an easy way to do it in R or in Python ?

Thank you in advance

Answer: Using **R**, *lag* function from *dplyr* package:

# example data
mybed <- read.table(text = "
1 3004104 3004110
1 3005819 3005825
1 3008315 3008321
1 3008893 3008899
1 3009812 3009818
1 3012422 3012428
1 3015794 3015800
1 3016183 3016189
1 3024019 3024025
1 3025279 3025285")

library(dplyr)

mybed %>%
transmute(
chr = V1,
start = lag(V2, default = 1),
end = V3)
Question: Hi,

Is it possible in plink to specify a non-human chromosome in the map file? I am curious if it can be run on a bacterial data set that only has 2 chromosomes and is haploid?

On the site it says:

By default, each line of the MAP file describes a single marker and must contain exactly 4 columns:

chromosome (1-22, X, Y or 0 if unplaced)
rs# or snp identifier
Genetic distance (morgans)
Base-pair position (bp units)

However, what if I just label everything chromosome 1 or 2?

Thanks!

Answer: If you read about the types of analysis that PLINK can do at the website, it's pretty clear that it's designed for humans. Or at least, diploid organisms which reproduce sexually, have chromosomes that undergo recombination and have sub-populations defined by SNPs. I'd also question whether bacteria even have SNPs in any meaningful sense.
Question: Hi.

I would like to search copy number variation using exome sequence data.

I have been using VarScan2 to detect copy number variations.

However, I am looking for more well-known or advanced variation in detecting copy number.

Can you suggest well-known and the most popular copy number variation tools within exome sequencing?

Answer: We use ExomeDepth ( http://cran.r-project.org/web/packages/ExomeDepth/index.html ) and CoNVeX ( http://www.biomedcentral.com/1471-2105/14/S2/S2 ) for copy number detection on targeted NGS panels (custom panel with bait density varied by depth of coverage from exome sequencing data), both have their pros and cons - with CoNVeX seemingly more susceptible to batch effects if there is alot of variation between samples.

We also use CoNVeX on whole exome data although you do get alot of false positives due to the high variablity in coverage.
Question: I want to get an evolutionary analysis by using NCBI NR database and Uniport database. But NR database has only blastdb and fasta format in which I can not find the taxon information. How can I get the NCBI NR database in genbank format.Thanks.

Answer: NCBI does not provide the NR database in GenBank format.

Alternatively, you could try the NCBI taxonomy database:

ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/

The taxdump.* are taxonomy database dumps and the detail introduction is available at: ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/taxdump_readme.txt.

The [accession2taxi][1]d folder provided the Accession number - taxonomy id relationships. You might also be interested in the gi_taxid_* files, which provide GI - taxonomy id relationship. But GenBank will not support GI soon.

[1]: ftp://ftp.ncbi.nlm.nih.gov/pub/taxonomy/accession2taxid/
Question: Can bedtools work directly on vcf.gz file, or do you have to first convert this into a vcf file?

Answer: Yes.

> bedtools will process gzipped BED, GFF, VCF and BEDPE files in the
> same manner as uncompressed files. Gzipped files are auto-detected
> thanks to a helpful contribution from Gordon Assaf.

http://bedtools.readthedocs.io/en/latest/content/overview.html
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: Cut and uniq to find duplicates, then grep -v them away. Relatively quick on 389MB dummy file.

time sort -V All_SNPs.txt | cut -f 1 | uniq -c | perl -ane 'if($F[0] ne "1"){print "$F[1]\t$F[0]\n";'} > All_SNPs.dup.chr-pos.txt

real 0m57.428s
user 3m19.091s
sys 0m3.152s

time cut -f 1 All_SNPs.dup.chr-pos.txt | grep -wvf - All_SNPs.txt > All_SNPs.nodup.txt

real 0m2.516s
user 0m1.697s
sys 0m0.229s
Question: I have RNA-seq data from the BALB/c mouse strain.

Looking for the reference genome on Ensembl, I found that, the most recent version, GRCm38 was build using the C57BL/6J strain.

I suppose that, the PATCH files contains haplotypes and variation also from other strain, like balb/c.

Are these information on the primary assembly file?

If I want to, How can I use the patch files for building the index? Just concatenate the files?

Or use the toplevel file for indexing?

Thank you.

Answer: If your goal is to reduce the allele bias in RNA-seq read mapping then using PATCH files won't be of much help. You will have to generate a strain-specific reference genome in order to align your reads or use some sensitive aligner. BALB/c has been sequenced as a part of Mouse Genome Project and the variant calls (VCF) can be downloaded from the following page: http://www.sanger.ac.uk/resources/mouse/genomes/. There are around 4 million SNPs and around 0.8 million indels between C57BL/6J and BALB/c. Although most of these variants fall into the intergenic regions but it would be a god practice to try to align reads in a haplotype-sensitive manner. You will have to create a customize reference genome by substituting small SNPs and Indels, and then perform the alignment. Sanger provides a big vcf file for all the strains, so you will have to 1) first extract variants for BALB/c strain and 2) then substitute them into the reference genome. There are many relevant posts on Biostars that have discussed both of these steps in a detailed manner.
Question: How should we select the correct kmer size for kmergenie analysis? I have read in wiki page that 'largest k-mer size to consider'. In that case I am looking for valuable comments from you people.

Answer: There may be a misunderstanding here.

"Selecting the correct kmer size for kmergenie analysis" implies that kmergenie would perhaps require user assistance?

This is precisely the opposite, it is designed to choose k for you. Just run kmergenie without any parameter, and it will explicitly print what it thinks is a good kmer size. Most of the time, the value is accurate.

Additionally, for unusual cases, you can take a look at the generated HTML report to make a more informed decision.
Question: I have a fastq file that seems to be contaminated by some sequences contaminating my reagents during library preparation. If I know the reads that came from reagents and I have them in a fasta format, do you think I can eliminate those reads from my fastq file? I want to remove any reads contaminating my fastq file. How can I work this out?

Answer: By using `bbduk.sh` [from BBMap][1]. Provide the contaminants as a multi-fasta file with `ref=` option.

[1]: http://seqanswers.com/forums/showthread.php?t=42776
Question: Hi!
Now I have 6 ATAC-seq samples which have been aligned to hg19 and called peaks by MACS2. 3 of them are from A cells and the other are from B cells. Now I'm going to do differential peak calling by MACS2. Before it, I need to merge peaks from 3 samples of the same cell type.
I use this command to call peaks.

macs2 callpeak -B -t sample_1.bam -n sample_1 --nomodel --extsize 120 --shift 60 \
--outdir sample_1 --keep-dup all

And then I get 5 output files.

sample_1_control_lambda.bdg
sample_1_peaks.narrowPeak
sample_1_peaks.xls
sample_1_summits.bed
sample_1_treat_pileup.bdg

But I'm confused how to merge them.
Shall I use *bedtools intersect* to merge 3 control_lambda.bdg files and treat_pileup.bdg files of A cell respectively? I'm worried if it will influence the value in the fourth column of the bdg files.

Or shall I merge their bam files before peak calling?

Answer: In my opinion "differential peak calling" does not exist with `macs2` as it is not able to make proper use of experimental replication and as such is not suitable to derive differential calls. It is great at calling peaks because it has been developed for that, but not more than that. I would call peaks on the merged files:

macs2 callpeak -t *.bam (...)

There is no need to merge the bam files as `t` can read multiple bams and then combine the pipeup internally.

Then convert the narrowPeak file to SAF format:

https://www.biostars.org/p/228636/

Then make a count matrix with `featureCounts` from the subread package:

featureCounts -a your.saf -F SAF -o out.counts *.bam

...and then feed this into something like DESeq2 or edgeR. Alternatively, look into the [DiffBind][1] or [csaw][2] packages which will instruct you on how to run a proper differential analysis.

[1]: https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf
[2]: http://bioconductor.org/books/3.13/csawBook/
Question: Hi,

I have a file like this:

> head(aa)
PLASER CURRELIG RTNPTHY
1 1 1 1
2 1 1 1
3 2 1 2
4 1 1 2
5 1 1 2
6 1 1 2

I have to create another column, say "new" which will satisfy these conditions:

PLASER=2 AND CURRELIG=1 -> new=2
RTNPHY=1 AND CURRELIG=1 -> new=1

I tried running this command:

aa$new <- ifelse(aa$CURRELIG==1 & aa$PLASER==2, 2, ifelse(aa$CURRELIG==1 & aa$RTNPHY==1,1,NA))

but I am getting this output which is not correct:

> head(aa)
PLASER CURRELIG RTNPTHY new
1 1 1 1 NA
2 1 1 1 NA
3 2 1 2 2
4 1 1 2 NA
5 1 1 2 NA
6 1 1 2 NA

Can someone please help with this?

Thanks
ANA

Answer: It looks like the `ifelse` command isn't meant to be nested like you have in your code, it's for binary output and it seem like you want three options (`1`, `2`, or `NA`). From the man page:

Usage:
ifelse(test, yes, no)
Arguments:
test: an object which can be coerced to logical mode.
yes: return values for true elements of ‘test’.
no: return values for false elements of ‘test’.

It might be easier to run:

aa$new <- NA
aa$new[aa$CURRELIG==1 & aa$PLASER==2] <- 2
aa$new[aa$CURRELIG==1 & aa$RTNPHY==1] <- 1

This will generate a new column of `NAs`, then replace the `NAs` with `1` when both `CURRELIG==1 & PLASER==2` or replace them with `2` when both `CURRELIG==1 & RTNPHY==1`.
Question: Hey I have a fasta file which contains more than 80 sequences. I want to take each sequences in different files and the file name should be header of corresponding sequence.

Answer: $ pip install pyfaidx
$ faidx --split-files in.fasta
Question: In terms of genotypes and strand information, would a sample genotype AG be equivalent to GA? And would genotype AG be equivalent to CT or TC for the same variant position?

Answer: Describing a variant on a different strand is always confusing.

If you are just talking about the genotype, you have to create the reverse complementary sequence. This means you start reading your initial genotype from right to left and changes each base to it complement (A/T or G/C).

So the reverse complement to AG would CT.

If you take the position of the variant into account, it looks a bit different. Positions are always described on the forward (+) strand. And if you describe an interval, you always starts with the lower position. So e.g. you have a variant AG on the forward strand starting on position 5, the equivalent description on the reverse strand would be TC starting on position 5 as well. You see I just created the complementary sequence and not the reverse complementary.

fin swimmer