Kezdőlap Felfedezés Súgó
Bejelentkezés
limour
/
limour-blog.github.io
tükrözi: https://dl.limour.top/cf/https://github.com/limour-blog/limour-blog.github.io.git
1
0
Másolás 0
Fájlok Problémák 0
Branch: main
Branch-ok Tag-ek
limour-blog.git...
/
source
/
_posts
/
2022-07-02-【迁移】STAR：一键脚本.md

2022-07-02-【迁移】STAR：一键脚本.md 14 KB
Állandó hivatkozás Előzmények Nyers

title: 【迁移】STAR：一键脚本 urlname: STAR--yi-jian-jiao-ben date: 2022-07-02 18:17:22 index_img: https://api.limour.top/randomImg?d=2022-07-02 18:17:22

tags: STAR

TCGA的STAR流程

{% note info %} 2023-09-27更新：添加了转录本计数相关内容 {% endnote %}

一些依赖

安装 biobambam

  • 项目地址：https://gitlab.com/german.tischler/biobambam2

    conda create -n biobambam -c bioconda biobambam -y

    安装 star

    conda create -n star -c bioconda star -y
conda activate star
conda install -c bioconda samtools -y
conda install -c bioconda rsem -y

可选：BAM转FASTQ

将所有bam文件以.bam结尾，单独存放到一个目录，conda activate biobambam，新建以下脚本运行

#!/bin/sh
ROOT=/home/jovyan/upload/22.07.02/muscle
mkdir $ROOT
 
for _ in *.bam
do
 
mkdir $ROOT/${_%.bam}
 
bamtofastq \
collate=1 \
exclude=QCFAIL,SECONDARY,SUPPLEMENTARY \
filename=$_ \
gz=1 \
inputformat=bam \
level=5 \
outputdir=$ROOT/${_%.bam} \
outputperreadgroup=1 \
outputperreadgroupsuffixF=_1.fq.gz \
outputperreadgroupsuffixF2=_2.fq.gz \
outputperreadgroupsuffixO=_o1.fq.gz \
outputperreadgroupsuffixO2=_o2.fq.gz \
outputperreadgroupsuffixS=_s.fq.gz \
tryoq=1 \
 
done
  • biobambam 高版本可能有 libmaus2 相关错误尚未修复
  • filename、outputdir 等参数等于号后不能有空格
  • 单端测序，outputdir 里只有 default_s.fq.gz 的输出
  • 更多信息

首次：构建小鼠的基因组索引

来源

  • human 和 mouse 推荐从 Gencode 上下载
  • 其他物种可以从 NCBI 上下载，也可以从 ENSEMBL 上下载
  • 需要参考基因组 fasta 文件和对应的 gtf 注释

步骤

zcat XL789vs123/XL1/default_s.fq.gz | head #这个例子只有50, 因此 sjdbOverhang 为49
 
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M29/GRCm39.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M29/gencode.vM29.primary_assembly.annotation.gtf.gz
 
gunzip *.gz
 
conda activate star
nano 2.sh
chmod +x 2.sh
 
#!/bin/sh
STAR \
--runMode genomeGenerate \
--genomeDir index \
--genomeFastaFiles GRCm39.primary_assembly.genome.fa \
--sjdbOverhang 49 \
--sjdbGTFfile gencode.vM29.primary_assembly.annotation.gtf \
--runThreadN 8
  • 如果要自己构建，可以使用 zcat R1.fq.gz | head 来查看reads长度，选用reads长度减1（149）作为 --sjdbOverhang 比默认的100要好，但是说明里认为绝大多数情况下100和理想值差不多

已构建好的索引

  • TCGA 上有已经构建好的索引，但是可能是损坏的
  • wget https://api.gdc.cancer.gov/data/07f2dca9-cd39-4cbf-90d2-c7a1b8df5139 -O star-2.7.5c_GRCh38.d1.vd1_gencode.v36.tgz
  • tar -zxvf star-2.7.5c_GRCh38.d1.vd1_gencode.v36.tgz

附加 构建转录本索引

cd /home/jovyan/upload/zl_liu/star/ # gtf文件和genome.fa文件所在目录
rsem-prepare-reference --gtf gencode.v36.primary_assembly.annotation.gtf GRCh38.primary_assembly.genome.fa GRCh38.primary_assembly.rsem.idx -p 8
grep -e '^ENST' GRCh38.primary_assembly.rsem.idx.ti > GRCh38.primary_assembly.rsem.t2n.ti
grep -e '^ENSG' GRCh38.primary_assembly.rsem.idx.ti > GRCh38.primary_assembly.rsem.g2n.ti

  • 大概生成了以下文件

    ├── [5.8M]  ./GRCh38.primary_assembly.rsem.g2n.ti
├── [ 770]  ./GRCh38.primary_assembly.rsem.idx.chrlist
├── [404K]  ./GRCh38.primary_assembly.rsem.idx.grp
├── [362M]  ./GRCh38.primary_assembly.rsem.idx.idx.fa
├── [362M]  ./GRCh38.primary_assembly.rsem.idx.n2g.idx.fa
├── [386M]  ./GRCh38.primary_assembly.rsem.idx.seq
├── [128M]  ./GRCh38.primary_assembly.rsem.idx.ti
├── [362M]  ./GRCh38.primary_assembly.rsem.idx.transcripts.fa
├── [6.5M]  ./GRCh38.primary_assembly.rsem.t2n.ti

    一键脚本

  • 原始数据质控

单独起一个目录，在这个目录下，以样本名为目录名，将样本的fataq文件以gzip压缩后存放，单端一个文件，双端两个文件，都可以，示例结构如下

XL1_12/
├── XL01
│   └── default_s.fq.gz
├── XL02
│   └── default_s.fq.gz
├── XL03
│   └── default_s.fq.gz
├── XL04
│   └── default_s.fq.gz
├── XL05
│   └── default_s.fq.gz
├── XL06
│   └── default_s.fq.gz
├── XL07
│   └── default_s.fq.gz
├── XL08
│   └── default_s.fq.gz
├── XL09
│   └── default_s.fq.gz
├── XL10
│   └── default_s.fq.gz
├── XL11
│   └── default_s.fq.gz
└── XL12
    └── default_s.fq.gz

在这个目录外，conda activate star，新建以下脚本运行，即可跑完前五步

#!/bin/bash
source activate star
#任务名
TASKN=XL1_12
#设置CleanData存放目录
CLEAN=/home/jovyan/upload/22.07.02/$TASKN
#设置输出目录
WORK=/home/jovyan/upload/22.07.02/output_$TASKN
#设置index目录
INDEX=/home/jovyan/upload/22.07.02/index
#设置参考文件位置
Reference=/home/jovyan/upload/22.07.02/GRCm39.primary_assembly.genome.fa
#设置 sjdbOverhang
sjdbOverhang=49
#设置 IIG 目录(这一步的输出目录)
IIG=/home/jovyan/upload/22.07.02/IIG_$TASKN
 
echo $CLEAN", "$WORK", "$INDEX
mkdir $WORK
 
CDIR=$(basename `pwd`)
echo $CDIR
echo $CLEAN
for file in $CLEAN/*
do
echo $file
SAMPLE=${file##*/}
echo $SAMPLE
mkdir $WORK"/"$SAMPLE
cd $WORK"/"$SAMPLE
 
STAR \
--genomeDir $INDEX \
--readFilesIn `ls $CLEAN/$SAMPLE/*` \
--runThreadN 4 \
--outFilterMultimapScoreRange 1 \
--outFilterMultimapNmax 20 \
--outFilterMismatchNmax 10 \
--alignIntronMax 500000 \
--alignMatesGapMax 1000000 \
--sjdbScore 2 \
--alignSJDBoverhangMin 1 \
--genomeLoad LoadAndRemove \
--readFilesCommand zcat \
--outFilterMatchNminOverLread 0.33 \
--outFilterScoreMinOverLread 0.33 \
--sjdbOverhang $sjdbOverhang \
--outSAMstrandField intronMotif \
--outSAMtype None \
--outSAMmode None \
 
done
 
echo $CLEAN", "$WORK", "$INDEX", "$IIG
mkdir $IIG
 
CDIR=$(basename `pwd`)
echo $CDIR
echo $CLEAN
for file in $CLEAN/*
do
echo $file
SAMPLE=${file##*/}
echo $WORK"/"$SAMPLE
done
 
STAR \
--runMode genomeGenerate \
--genomeDir $IIG \
--genomeFastaFiles $Reference \
--sjdbOverhang $sjdbOverhang \
--runThreadN 4 \
--sjdbFileChrStartEnd `ls $WORK/*/SJ.out.tab`
 
ln -s $INDEX/exonGeTrInfo.tab $IIG
ln -s $INDEX/exonInfo.tab $IIG
ln -s $INDEX/geneInfo.tab $IIG
ln -s $INDEX/sjdbList.fromGTF.out.tab $IIG
ln -s $INDEX/transcriptInfo.tab $IIG
 
CDIR=$(basename `pwd`)
echo $CDIR
echo $CLEAN
for file in $CLEAN/*
do
echo $file
SAMPLE=${file##*/}
echo $WORK"/"$SAMPLE
mkdir $WORK"/"$SAMPLE"/Res"
cd $WORK"/"$SAMPLE"/Res"
 
STAR \
--genomeDir $IIG \
--readFilesIn `ls $CLEAN/$SAMPLE/*` \
--runThreadN 8 \
--quantMode TranscriptomeSAM GeneCounts \
--outFilterMultimapScoreRange 1 \
--outFilterMultimapNmax 20 \
--outFilterMismatchNmax 10 \
--alignIntronMax 500000 \
--alignMatesGapMax 1000000 \
--sjdbScore 2 \
--alignSJDBoverhangMin 1 \
--genomeLoad LoadAndRemove \
--limitBAMsortRAM  35000000000 \
--readFilesCommand zcat \
--outFilterMatchNminOverLread 0.33 \
--outFilterScoreMinOverLread 0.33 \
--sjdbOverhang $sjdbOverhang \
--outSAMstrandField intronMotif \
--outSAMattributes NH HI NM MD AS XS \
--outSAMunmapped Within \
--outSAMtype BAM SortedByCoordinate \
--outSAMheaderHD @HD VN:1.4 \
--outSAMattrRGline ID:sample SM:sample PL:ILLUMINA
 
done

附加 转录本计数

  • 转录本水平: cufflinks、eXpress、RSEM

  • 外显子水平: DEXseq

    nano RSEM.sh && chmod +x RSEM.sh
./RSEM.sh

    #!/bin/bash
source activate star
 
#任务名
TASKN=cleandata
#设置上一步的输出目录
WORK=/home/jovyan/work/shRNA-seq/output_$TASKN
#设置参考文件位置
Reference=/home/jovyan/upload/zl_liu/star/GRCh38.primary_assembly.genome.fa
#设置转录本索引位置(前缀)
RSEM_IDX=/home/jovyan/upload/zl_liu/star/GRCh38.primary_assembly.rsem.idx
#设置这一步的输出目录
RSEM_OUT=/home/jovyan/work/shRNA-seq/RSEM_$TASKN
mkdir -p $RSEM_OUT
 
for file in $WORK/*
do
echo $file
SAMPLE=${file##*/}
echo $SAMPLE
mkdir $RSEM_OUT"/"$SAMPLE
cd $RSEM_OUT"/"$SAMPLE
 
rsem-calculate-expression -no-bam-output --alignments -p 8 --paired-end \
$WORK"/"$SAMPLE"/Res/Aligned.toTranscriptome.out.bam" \
$RSEM_IDX \
'resm'
 
done

  • 最后得到的RSEM_OUT目录结构如下

    [158M]  .
├── [ 26M]  ./NC-1
│   ├── [6.8M]  ./NC-1/resm.genes.results
│   ├── [ 15M]  ./NC-1/resm.isoforms.results
│   └── [4.3M]  ./NC-1/resm.stat
│       ├── [ 871]  ./NC-1/resm.stat/resm.cnt
│       ├── [493K]  ./NC-1/resm.stat/resm.model
│       └── [3.9M]  ./NC-1/resm.stat/resm.theta
├── [ 26M]  ./NC-2
│   ├── [6.8M]  ./NC-2/resm.genes.results
│   ├── [ 15M]  ./NC-2/resm.isoforms.results
│   └── [4.4M]  ./NC-2/resm.stat
│       ├── [ 872]  ./NC-2/resm.stat/resm.cnt
│       ├── [493K]  ./NC-2/resm.stat/resm.model
│       └── [3.9M]  ./NC-2/resm.stat/resm.theta
├── [ 26M]  ./NC-3
│   ├── [6.8M]  ./NC-3/resm.genes.results
│   ├── [ 15M]  ./NC-3/resm.isoforms.results
│   └── [4.4M]  ./NC-3/resm.stat
│       ├── [ 872]  ./NC-3/resm.stat/resm.cnt
│       ├── [493K]  ./NC-3/resm.stat/resm.model
│       └── [3.9M]  ./NC-3/resm.stat/resm.theta
├── [ 26M]  ./shRNA1-1
│   ├── [6.8M]  ./shRNA1-1/resm.genes.results
│   ├── [ 15M]  ./shRNA1-1/resm.isoforms.results
│   └── [4.4M]  ./shRNA1-1/resm.stat
│       ├── [ 874]  ./shRNA1-1/resm.stat/resm.cnt
│       ├── [493K]  ./shRNA1-1/resm.stat/resm.model
│       └── [3.9M]  ./shRNA1-1/resm.stat/resm.theta
├── [ 26M]  ./shRNA1-2
│   ├── [6.8M]  ./shRNA1-2/resm.genes.results
│   ├── [ 15M]  ./shRNA1-2/resm.isoforms.results
│   └── [4.3M]  ./shRNA1-2/resm.stat
│       ├── [ 849]  ./shRNA1-2/resm.stat/resm.cnt
│       ├── [493K]  ./shRNA1-2/resm.stat/resm.model
│       └── [3.8M]  ./shRNA1-2/resm.stat/resm.theta
└── [ 26M]  ./shRNA1-3
├── [6.8M]  ./shRNA1-3/resm.genes.results
├── [ 15M]  ./shRNA1-3/resm.isoforms.results
└── [4.3M]  ./shRNA1-3/resm.stat
    ├── [ 867]  ./shRNA1-3/resm.stat/resm.cnt
    ├── [493K]  ./shRNA1-3/resm.stat/resm.model
    └── [3.8M]  ./shRNA1-3/resm.stat/resm.theta

  • 组装RSEM_Counts文件

    #设置转录本索引转换的位置(前缀)(无idx)
RSEM_x2n='/home/jovyan/upload/zl_liu/star/GRCh38.primary_assembly.rsem'
#设置RSEM的输出目录
RSEM_OUT='/home/jovyan/work/shRNA-seq/RSEM_cleandata'
file_list <- list.dirs(RSEM_OUT, recursive=F, full.names = F)
t2n <- read.table(file = paste0(RSEM_x2n, '.t2n.ti') , sep = '\t', header = F)
g2n <- read.table(file = paste0(RSEM_x2n, '.g2n.ti') , sep = '\t', header = F)
transcript_counts <- cbind(t2n, g2n)
colnames(transcript_counts) <- c('transcript_id', 'transcript_name', 'gene_id', 'symbol')
transcript_tpm <- transcript_counts
transcript_IsoPct <- transcript_counts
for(sample in file_list){
test_tab <- read.table(file = file.path(RSEM_OUT, sample, 'resm.isoforms.results') , sep = '\t', header = T)
transcript_counts[sample] <- test_tab$expected_count
transcript_tpm[sample] <- test_tab$TPM
transcript_IsoPct[sample] <- test_tab$IsoPct
}
write.csv(x = transcript_counts, file = 'transcript_counts.csv') # 期望计数
write.csv(x = transcript_tpm, file = 'transcript_tpm.csv')
write.csv(x = transcript_IsoPct, file = 'transcript_IsoPct.csv') # 该转录本占所属基因的比例

    组装的转录本计数文件格式

    组装Counts文件

    #设置第三步的输出目录
WORK='/home/jovyan/upload/22.07.02/output_XL789vs123'
#设置index中基因注释位置
Reference='/home/jovyan/upload/22.07.02/index/geneInfo.tab'
file_list <- list.dirs(WORK, recursive=F, full.names = F)
geneN <- read.table(file = Reference, sep = '\t', skip = 1)
colnames(geneN) <- c('ID', 'symbol', 'type')
for(sample in file_list){
test_tab <- read.table(file = file.path(WORK, sample, 'Res', 'ReadsPerGene.out.tab') , sep = '\t', header = F)
test_tab <- test_tab[-c(1:4), ]
geneN[sample] <- test_tab[2]
}
write.csv(x = geneN, file = 'XL789vs123.csv')

组装的Counts文件格式

DESeq2分析

library(DESeq2)
count_all <- read.csv("liver.csv",header=TRUE,row.names=1)
count_all
cts_b <- count_all[ ,c(-1,-2,-3)]
rownames(cts_b) <- count_all$ID
cts_bb <- cts_b
 
cts_b <- cts_bb[,c('XL07', 'XL08', 'XL09', 'XL10', 'XL11', 'XL12')]
keep <- rowSums(cts_b) > 10
cts_b[keep,]
conditions <- factor(c(rep("Control",3), rep("XL",3)))    
colData_b <- data.frame(row.names = colnames(cts_b), conditions)
colData_b
dds <- DESeqDataSetFromMatrix(countData = cts_b[keep,],
                              colData = colData_b,
                              design = ~ conditions)
dds <- DESeq(dds)
res <- results(dds)
rres <- cbind(count_all[keep,c(1,2,3)], data.frame(res))
write.csv(rres, file='XL101112vs789_DESeq2.csv')
rres

附加：删除临时文件

#!/bin/sh
#设置CleanData存放目录
CLEAN=/home/jovyan/upload/yy_zhang（备份）/RNA-seq/Cleandata
#设置这一步的输出目录 (确保目录存在)
WORK=/home/jovyan/upload/zl_liu/star_data/yyz_01/output
 
for file in $CLEAN/*
do
echo $file
SAMPLE=${file##*/}
echo $SAMPLE
 
rm -rf $WORK"/"$SAMPLE"/IIG"
 
done