Salmon代码备存

#Salmon代码备存

听说你还可以分析scRNA

salmon alevin -l ISR -1 Pi_F_S1_L005_R1_001.fastq.gz -2 Pi_F_S1_L005_R2_001.fastq.gz --chromium -i ~/index/ggzjn_index/ -p 10 --dumpMtx -o alevin_output/ --tgMap txp2gene.txt --expectCells 18000
salmon alevin -l ISR -1 Pi_M_S1_L004_R1_001.fastq.gz -2 Pi_M_S1_L004_R2_001.fastq.gz --chromium -i ~/index/ggzjn_index/ -p 15 --dumpMtx -o alevin_output/ --tgMap txp2gene.txt --expectCells 18000

salmon alevin -l ISR -1 LI_F_v3_S1_L001_R1_001.fastq.gz -2 LI_F_v3_S1_L001_R2_001.fastq.gz --chromiumV3 -i ~/index/ggzjn_index/ -p 10 --dumpMtx --dumpFeatures -o LI_F_v3_alevin_output/ --tgMap txp2gene.txt --expectCells 25000

salmon alevin -l ISR -1 LI_M_v3_S1_L001_R1_001.fastq.gz -2 LI_M_v3_S1_L001_R2_001.fastq.gz --chromiumV3 -i ~/index/ggzjn_index/ -p 10 --dumpMtx --dumpFeatures -o LI_M_v3_alevin_output/ --tgMap txp2gene.txt --expectCells 25000

数据下载-plan A

#!/bin/bash
for fn in SRR1299{036..049};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
ascp -QT -l 300m -P33001 -i ~/miniconda3/etc/asperaweb_id_dsa.openssh [email protected]:/vol1/srr/SRR707/008/${fn} ./
mv ${fn} ${fn}.sra
done

数据下载-plan B

nano srr.txt

SRR1524238
SRR1524239
SRR1524240
SRR1524241

prefetch --option-file srr.txt --max-size 100G -O ./
prefetch -t ascp -a "/home/labwang/.aspera/connect/bin/ascp|/home/labwang/.aspera/connect/etc/asperaweb_id_dsa.openssh" --option-file srr.txt -O ./
# 啊，升级了好麻烦啊！！！ 
prefetch -t ascp -a "/home/labwang/.aspera/connect/bin/ascp|/home/labwang/.aspera/connect/etc/asperaweb_id_dsa.openssh" --option-file srr.txt
cd /usr/local/ncbi/sra-tools/bin  然后 ./vdb-config -i 然后 巴拉巴拉~~~

分析

salmon index -t Taeniopygia_guttata.taeGut3.2.4.cdna.all.fa -i finch2_index

#!/bin/bash
for fn in SRR7072{158..175};
do
path=$(pwd)
samp=`basename ${fn}`
echo "Processing sample ${samp}"
pfastq-dump -t 20 --split-files --gzip -O ${path}/fastq ${samp}.sra
done

#!/bin/bash
for fn in SRR7072{158..175};
do
path=$(pwd)
echo "Processing sample ${fn}"
salmon quant -i ~/index/ggzjn_index --validateMappings -l A \
         -1 ${path}/fastq/${fn}_1.fastq.gz \
         -2 ${path}/fastq/${fn}_2.fastq.gz \
         -p 20 -o ${path}/quant/${fn}_quant
done

#!/bin/bash
for filename in *_*;
do
path=$(pwd)
samp=`basename ${filename}`
echo "Processing sample ${filename}"
salmon quant -i ~/index/ggzjn_index --validateMappings -l A \
        -1 ${path}/${samp}/${samp}_1.fq.gz \
        -2 ${path}/${samp}/${samp}_2.fq.gz \
        -p 10 -o ${path}/quant/${samp}_quant
done

bash /data3/PI-LV/R/DESEQ2.sh -f result-DE.txt -s samplefiles.txt -p compare_pair.txt 

bash /data3/PI-LV/R/DESEQ2.sh -f result-DE.txt -s samplefiles.txt -p compare_pair.txt -F 0.5 -P 0.01

library("tximport")
library(readr)

dir <- "D:/Users/zjn-lab/Documents/salmon"
list.files(dir)
sample <- paste0("ERR348",seq(560,588),"_quant")
files <- file.path(dir,"quant",sample,"quant.sf")
names(files) <- paste0("sample",c(1:3))
all(file.exists(files))
head(files)

tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v97.csv"))
head(tx2gene)


txi <- tximport(files, type = "salmon", tx2gene = tx2gene)
head(txi)
names(txi)
head(txi$abundance)
head(txi$length)
head(txi$counts)
head(txi$countsFromAbundance)

#write.table(txi$abundance, "abundance.txt", sep="\t", quote=F)
write.csv(txi$abundance, "abundance.csv", quote=F)


txi.tx <- tximport(files, type = "salmon", txOut = TRUE)
txi.sum <- summarizeToGene(txi.tx, tx2gene)
all.equal(txi$counts, txi.sum$counts)
head(txi.tx)
names(txi.tx)
head(txi.tx$abundance)
head(txi.tx$length)
head(txi.tx$counts)
head(txi.tx$countsFromAbundance)

library(DESeq2)

sampleTable <- data.frame(condition = factor(rep(c("untreated", "treated"), each = 4)))
rownames(sampleTable) <- colnames(txi$counts)
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~condition)
dds
sampleTable

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,] 

数据质检过滤

#!/bin/bash
for fn in SRR1524{238..249};
do
path=$(pwd)
samp=`basename ${fn}`
echo "Processing sample ${samp}"
fastqc -o ${path}/fastqc -t 50 ${path}/fastq/${fn}_1.fastq.gz
fastqc -o ${path}/fastqc -t 50 ${path}/fastq/${fn}_2.fastq.gz
done

multiqc .

#!/bin/bash
for fn in SRR1524{238..249};
do
samp=`basename ${fn}`
echo "Processing sample ${samp}"
trimmomatic PE -phred33 -threads 20 ${fn}_1.fastq.gz ${fn}_2.fastq.gz ${fn}_paired_1.fastq.gz ${fn}_unpaired_1.fastq.gz ${fn}_paired_2.fastq.gz ${fn}_unpaired_2.fastq.gz ILLUMINACLIP:/home/labwang/miniconda3/share/trimmomatic-0.38-1/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
rm -r ${fn}_unpaired_1.fastq.gz
rm -r ${fn}_unpaired_2.fastq.gz
rm -r ${fn}_1.fastq.gz
rm -r ${fn}_2.fastq.gz
done

去头去尾

trimmomatic PE -phred33 -threads 20 ${fn}_1.fq.gz ${fn}_2.fq.gz ${fn}_paired_1.fastq.gz ${fn}_unpaired_1.fastq.gz ${fn}_paired_2.fastq.gz ${fn}_unpaired_2.fastq.gz CROP:89 HEADCROP:10 MINLEN:36

update

#!/bin/bash
# get all filename in specified path
path=$(pwd)
files=$(ls $path)
for filename in $files
do
echo $filename
path=$(pwd)
salmon quant -i ~/index/ggzjn_index --validateMappings -l A \
-1 ${path}/${filename}/${filename}_1.fq.gz \
-2 ${path}/${filename}/${filename}_2.fq.gz \
-p 10 -o /data3/lv-hy/quant/${filename}_quant
done