search

ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
124 stars 31 forks source link

不太一样的Bulk转录组上游分析 #2515

Closed ixxmu closed 2 years ago

ixxmu commented 2 years ago

https://mp.weixin.qq.com/s/dz3ImtuGdZEs27TTBda6VQ

github-actions[bot] commented 2 years ago

不太一样的Bulk转录组上游分析 by 生信菜鸟团

测序类型：猪，转录组双端测序数据

问题：我们自己鉴定到一个全新的基因，人为命名为Novel_Genes，测序获得了它的序列（fa文件），但是NCBI和ENSEMBL等数据库并没有收录这个基因ID。如何从RNA-seq中获取其表达量？

分析思路：常规的转录组上游分析得到bam文件；然后将该Novel_Genes的fa转化为gtf文件，进行featureCount分析。

  • 用于示例的目标基因的序列（这里的序列是我随机挑选用于举例的，并无实际意义），我将其命名为Novel_Genes：
cat test.fa

Step1. 环境部署及fastq文件下载

参考 Linux随笔（二）“软件管家“ conda

  • conda和mamba下载

  • 环境配制

    # 建立conda环境
    conda env create -n rnaseq python==3
    # 激活conda环境
    source activate rnaseq
    # 安装软件 
    mamba install -y fastqc multiqc hisat2 samtools subread  trim-galore
    mamba install -y -c bioconda blast #用于blast

  • 双端fastq测序数据下载：这里已经下载好了。fastq下载我们可以使用aspera加速下载，参考 91例免疫治疗队列（含生存）转录组上游分析实战

  • 猪的基因组文件下载自Ensembl数据库：

#参考基因组
wget -c http://ftp.ensembl.org/pub/release-107/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz
gunzip Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz
#注释文件
wget -c http://ftp.ensembl.org/pub/release-107/gtf/sus_scrofa/Sus_scrofa.Sscrofa11.1.107.chr.gtf.gz
gunzip Sus_scrofa.Sscrofa11.1.107.chr.gtf.gz
  • 构建bowtie2索引文件：
mkdir index
cd index
# bowtie2构建索引
bowtie2-build --threads 10 ../Sus_scrofa.Sscrofa11.1.dna.toplevel.fa Sus_scrofa

会生成6个bt2文件

ls -lh *bt2

Step2.Blast check

文件夹部署：

mkdir 1.fastq_data  2.blast 3.fastqc 4.clean 5.mapping 6.count

先Blast一下这个目的基因的序列，参考https://blog.csdn.net/hac_kill_you/article/details/121616212

cd 2.blast
# 将pig.fastq作为reference建立库，然后用test.fa比对
# 建立库 （将其中一条序列建库，作为reference）
makeblastdb -in ../Sus_scrofa.Sscrofa11.1.dna.toplevel.fa -dbtype nucl -out pig_db

# 比对 （test.fa作为query）
# -db 就是上一步建库的名字（-out的参数）
blastn -query ../test.fa -db pig_db -out result

查看blast的结果：

cat result
BLASTN 2.12.0+

Reference: Zheng Zhang, Scott Schwartz, Lukas Wagner, and Webb
Miller (2000), "A greedy algorithm for aligning DNA sequences", J
Comput Biol 2000; 7(1-2):203-14.

Database: ../Sus_scrofa.Sscrofa11.1.dna.toplevel.fa
           613 sequences; 2,501,912,388 total letters

Query= Novel_Genes

Length=809
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

2 dna:chromosome chromosome:Sscrofa11.1:2:1:151935994:1 REF           1495    0.0  

>2 dna:chromosome chromosome:Sscrofa11.1:2:1:151935994:1 REF
Length=151935994

 Score = 1495 bits (809),  Expect = 0.0
 Identities = 809/809 (100%), Gaps = 0/809 (0%)
 Strand=Plus/Minus

Query  1        CGCTCACTGACCCGCATCTGACACCTGATACCCAGCGGCCAGCAGAGGTGTGGCCCCATA  60
                ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  6114196  CGCTCACTGACCCGCATCTGACACCTGATACCCAGCGGCCAGCAGAGGTGTGGCCCCATA  6114137

......

Query  781      AATGGCCGCAGAGTCCAGGCCAGCGTGAG  809
                |||||||||||||||||||||||||||||
Sbjct  6113416  AATGGCCGCAGAGTCCAGGCCAGCGTGAG  6113388

由Blast的结果可知，此Novel_Genes在chr2:6113388-6114196区间，为反义链，全长809bp

制作此Novel_Genes的GTF文件

echo -e '2\tensembl\texon\6113388\6114196\t.\t-\t.\tgene_id "Novel_Genes"; transcript_id "Novel_Genes"; gene_name "Novel_Genes"; gene_biotype "lncRNA";' >Novel_Genes.gtf

将此gtf添加到Sus_scrofa.Sscrofa11.1.107.chr.gtf

#做一个备份
cp Sus_scrofa.Sscrofa11.1.107.chr.gtf self.add.Novel_Genes.gtf
#追加
cat Novel_Genes.gtf >>self.add.Novel_Genes.gtf
#check
cat self.add.Novel_Genes.gtf | tail

Step3.常规的转录组上游分析

3.1 fastqQC质控

cd 3.fastqc
# 使用FastQC软件对单个fastq文件进行质量评估
ln -s ../1.fastq_data/*.fastq.gz ./

###生成的fastqc放到3.fastqc中,-t指定线程数 后台运行
nohup fastqc -t 5 -o ./ ./*.fastq.gz &

###使用MutliQC整合FastQC结果，将后缀为.html的文件进行multiqc处理
multiqc ./

3.2 clean: trim_galore进行过滤

cd ../4.clean/

###做一个软连接文件,将格式转换后的fq.gz文件链接至4.clean
ln -s ../1.fastq/*.fastq.gz ./

##########i=${i/_1.fastq.gz/}意思是除去“i”后面的“_1.fastq.gz”
for i in `ls *_1.fastq.gz`
do
i=${i/_1.fastq.gz/}
echo "trim_galore --phred33 -q 25 --cores 8 --length 36 -e 0.1 --stringency 3 --paired -o ./ ${i}_1.fastq.gz ${i}_2.fastq.gz"
done > trim_galore.command

# multiqc
nohup bash trim_galore.command 1>trim_galore.log 2>&1 &

#批量kill
#ps -ef | grep trim_galore | awk '{print $2}' | while read id;do kill $id;done  

fastq_results质控：

###生成的fastqc放到~/sra/3.fastq_qc/中,-t指定线程数。
nohup fastqc -t 30 -o ./ ./*_trimmed.fq.gz & #单个单个文件进行操作

###使用MutliQC整合FastQC结果，将后缀为.html的文件进行multiqc处理
multiqc ./*.zip

3.3 比对

这里是猪的双端测序数据，index下载及构建见开头部分。

cd ../5.mapping
ln -s ../4.clean/*.fastq.gz ./

## 比对到猪参考基因组,并用samtools的sort排一下序
cat > bowtie2_paried.command
index=~/tmp_data/pig_fastq/index/Sus_scrofa
for id in `ls *_1.fastq.gz`
do
file=$(basename $id)
sample=${file%%_1.fastq.gz}
bowtie2 -p 8 -x ${index} -1 ${sample}_1.fastq.gz -2 ${sample}_2.fastq.gz | samtools sort -O bam -@ 8 -o - >${sample}.bam
done 

###运行脚本
nohup bash bowtie2_paried.command &
# 查看任务进度
ps -ef |grep bowtie2
#ps -ef | grep bowtie2 | awk '{print $2}' | while read id;do kill $id;done  #批量Kill

期间有一个样本报错了：

Error: Read SRR429934.9367055 HIGAATCGTAGGTGG044/1 has more quality values than read characters.
(ERR): bowtie2-align died with signal 6 (ABRT) (core dumped)
[bam_sort_core] merging from 0 files and 10 in-memory blocks...

我怀疑是样本没下载完整，我们check一下：

function check_gzip_file {
 local file=$1
 gzip -t $file && echo ok || echo bad
}

check_gzip_file SRR4299265_1.fastq.gz
check_gzip_file SRR4299265_2.fastq.gz

#其实一开始就应该check所有的文件
for id in `ls *.fastq.gz`
do
file=$(basename $id)
echo $file
check_gzip_file $file 
done > check.log

果然是有一个文件没有下载完整。这里只要重新下载跑上述流程即可:

ascp -v -QT -l 300m -P33001 -k1 \
-i ~/miniconda3/envs/download/etc/asperaweb_id_dsa.openssh   \
era-fasp@fasp.sra.ebi.ac.uk:/vol1/fastq/SRR429/005/SRR4299265/SRR4299265_1.fastq.gz .

所有fastq文件转为bam文件后，构建index并质控查看比对成功率：

## 构建index
ls *bam |xargs -i samtools index {} -@ 10

## QC
ls *bam |while read id ;do (nohup samtools flagstat -@ 1 $id > $(basename $id ".bam").flagstat &);done

##查看比对成功率
grep "properly paired" *.flagstat

Step.4 featureCount定量

计数的软件有很多：htseq / featureCounts / bedtools / Deeptools / salmon

featureCount是subread套件的一个模块，最大的优点就是速度非常快，使用全部overlap的reads计数，灵活考虑多比对的reads的计数。

这里使用上述自定义构建的猪的gtf注释（含目的Novel_Genes）：self.add.Novel_Genes.gtf

cd ../6.count
### 使用sort后的bam文件进行操作
ln -s ../5.mapping/*.bam ./

### reads计数，其中的-t,-g都是默认的。其中-g可以指定显示gtf文件中attributes那一列中的任意值。（例：gene_id，gene_name，transcript_id等）
### -T 45 四十五个线程；-p 双端测序；-t exon是要统计的feature的名称；-g gene_id是要统计的meta-feature的名称
nohup featureCounts -T 45 -p -t exon -g gene_id -a ../self.add.Novel_Genes.gtf  -o ./all.id.txt *.bam &

### 简化结果，去除特征列，只保留基因计数的列
cat all.id.txt | cut -f1,7-1000 > counts.txt

# 查看任务进度
ps -ef | grep featureCounts
ps -ef | grep featureCounts | awk '{print $2}' | while read id;do kill $id;done  #批量Kill

可以看到自定义的这个Novel_Genes基因的表达量在最后一列，其他标准基因都是ENSEMBL ID：

cat counts.txt | tail

一个小技巧，希望能够抛砖引玉～

- END -