6

Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose. C++, Rust, Go, Nim and Julia can be as fast. // Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h // Compile with: gcc -O2 -o myprog this-file.c -lz #include <zlib.h> ...


6

To answer your direct question, there are a few reasons why there might be high levels of sequence duplication. From the FastQC help: The underlying assumption of this module is of a diverse unenriched library. Any deviation from this assumption will naturally generate duplicates and can lead to warnings or errors from this module. As @DevonRyan ...


6

FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether ...


5

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools. But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as ...


4

With bwa-mem or minimap2, the recommended way is samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p - This passes the BC and RX tags through bwa-mem and copies them to the SAM output. If you are running paired-end alignment and your BAM is coordinate sorted, you have to run samtools collate to group reads by name: samtools collate -Ol0 in.bam tmp | ...


4

I wrote a python script based on pysam that is free for anyone to use: """ MIT License Copyright (c) 2020 Warren W. Kretzschmar Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights ...


4

EDIT after OP's update: Try this Perl one-liner: perl -pe 'BEGIN { $i = 1 } $i += s{>([^_]+)_.*_}{>RNA${i}#${1}/}' input_file > output_file Here, the command line flags are: -e: tells the perl interpreter to look for the code inline, rather than in a file with the script, -p: loop over the input one line at a time, execute the code in the one-...


3

Your own answer doesn't seem to do what you ask as it doesn't add the sequential identifier to each sequence. Maybe try this? awk -v rec=0 '{ if($0 ~ "^>"){ rec++ sub(">", ">RNA" rec "#", $0) sub("_.*_", "/", $0); } print $0 }' yourfile.fa


3

In my experience, it is more common to use all clinical data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data. If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.


3

From the link you posted it looks as though they used UMI tools as per the section which says 'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'. The link above has instructions for installation and use. From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read. Looking ...


3

The command where you trim with adapters and by quality is perfectly fine. That FastQC isn't perfectly happy is expected. It's a tool made for whole genome sequencing QC, you SHOULD see a number of "fail"s with any RNA-seq protocol. Steps that should always fail in RNA-seq: per-base sequence content (there's "random priming" performed that isn't completely ...


3

In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the ...


3

I don't know if this question has been solved already, but what they try to do is equalize the depth of sequencing for each cell. Therefore, they scale for the total number of reads. If you regress out (via linear or negative binomial regression) the differences in the number of reads per cell, you end up with cells that have been sequenced with the same ...


2

Based on the info you provide, ArrayBin R package provides you the necessary tools: binarize.array() from ArrayBin, allowing: Implementation of an adaptive method for binarizing gene expression data on a per-probe basis and demonstrate the superior effectiveness of our method when compared with other, commonly used approaches. This adaptive binarization ...


2

I am not sure how much you know about bioinformatics already, can you use R? For a bioinformatician looking at QC for microarrays should not be a big deal, at least for me it would take maybe a day (or two) to get this done. However, if you never used R and want to start from scratch, it depends on how quickly you learn how to deal with arrays and QC. There ...


2

You can filter records in Python for example using pyvcf.VcfFrame.filter_flagany() method I wrote: >>> from fuc import pyvcf >>> data = { ... 'CHROM': ['chr1', 'chr1', 'chr1', 'chr1'], ... 'POS': [100, 101, 102, 103], ... 'ID': ['.', '.', '.', '.'], ... 'REF': ['G', 'T', 'A', 'C'], ... 'ALT': ['A', 'C', 'T', 'A'], ... ...


2

sed -e 's:^\(.*\)_.*_\(.*\)$:\1/\2:' < input > output


2

In A2M format, upper case letters represent matches, lower case letters represent inserts, dashes represent deletions, and dots (or spaces outside the identifier lines) represent gaps aligned to inserts. So, both "-" and "." are essentially gaps, but assumed to have different origin. This information is supplementary, and most MSA ...


1

It was done using sangerseqR with rpy2 in python import rpy2.robjects as r from rpy2.robjects.packages import importr utils = importr('utils') utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR") utils.chooseBioCmirror(ind=1) # select the first mirror in the list utils....


1

Depending on the cluster management tool, you might have received e-mails when the "job" begins and ends. If so, you can check the "Exit status" of the job. For example, in the case of our HPC the relevant lines from the "job completion e-mail" would be: Execution terminated Exit_status=0 BTW, at first I was confused with the phrase "Execution terminated", ...


1

Removing outliers is common practice in statistical modeling and perfectly acceptable. However, with regards 1.5 IQR I am far from certain about this approach. Normally, if you want to be conservative then 3 standard deviations (SD) denote an outlier, which is more stringent than IQR. Some use 2 SD. If the value is lower than 2 SD from the group mean it isn'...


1

Hard to say without seeing VCF rows, but generally speaking with text files you can use grep, if I understand the question: grep -vw "Match=EXACT" file.vcf > filtered.vcf


1

Assuming your "master" files looks like this: df <- data.frame(Name = LETTERS[1:5], Col1 = sample(1:10,5), Col2 = sample(1:10,5), Col11 = sample(1:10,5), Col3 = sample(1:10,5)) Name Col1 Col2 Col11 Col3 1 A 5 4 3 10 2 B 1 6 5 6 3 C 3 1 ...


1

Example in R. If you have a file of column names, names.txt, such as: mpg cyl disp and a data.frame like mtcars > mtcars mpg cyl disp hp drat wt qsec vs am gear carb Mazda RX4 21.0 6 160.0 110 3.90 2.620 16.46 0 1 4 4 Mazda RX4 Wag 21.0 6 160.0 110 3.90 2.875 17.02 0 1 4 4 Datsun 710 22.8 ...


1

There's a million different ways to do simple taks like this. Pick one and learn it:) Python is my go to as it can do pretty much everything you'd ever want. I saved your example as csv with no spaces called expression.txt. import pandas as pd df = pd.read_csv('expression.txt') df[['SM-5GZZ7','SM-5GIEN','SM-5EGGH']] SM-5GZZ7 SM-5GIEN SM-5EGGH 0 0 ...


1

In general you want to: Clean your text Split the features by semicolons Split keys and values by '=' Create a temporary data frame with keys and value Reshape library(tidyverse) library(urltools) library(textclean) test_df <- data.frame(pos=c(3002983, 3002881, 3010946, 3021775), metadata=c("ID=SRX661585;Name=Ctcf%20(@%20Brain);Title=GSM1446329:%...


1

Via BEDOPS convert2bed (psl2bed) and bedops operations: $ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed The file answer.bed will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases. This BED file can be run through samtools faidx or similar to get sequence data. References: https://bedops....


1

I wrote a package called sinto (https://github.com/timoast/sinto) that contains a function to split a BAM file into different files based on the cell barcodes. See the documentation here: https://timoast.github.io/sinto/basic_usage.html#filter-cell-barcodes-from-bam-file


Only top voted, non community-wiki answers of a minimum length are eligible