【转录组04】参考基因组

使用两个软件对fq数据进行比对，得到比对文件sam/bam，并探索比对结果。

• 三个常用参考基因组数据库：
  Ensembl：www.ensembl.org
  NCBI：
  https://www.ncbi.nlm.nih.gov/projects/genome/gu
  ide/human/index.shtml
  UCSC：http://www.genome.ucsc.edu/

• Ensembl官网

  
  image.png

image.png
image.png 这种ftp类型的地址是可以直接用weget命令进行下载的 image.png

## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：http://asia.ensembl.org/index.html
# ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/

# 进入到参考基因组目录
cd /teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95

# 下载基因组序列
wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

——————————————————数据太大就不下载了——————————使用现成的————————————————————————————-

软链接下载好的参考基因组到自己的目录下

• 参考基因组注释文件

  
  每一行有九列

image.png image.png image.png

# 下载基因组注释文件
wget ftp://ftp.ensembl.org/pub/release-95/gtf/homo_sapiens/Homo_sapiens.GRCh38.95.gtf.gz

习题

1.fastq与fasta文件转换

应用：加深对两种文件格式的理解，并且fa有后续应用。

#进入fastq目录
less -S SRR1039510_1.fastq.gz | paste - - - - |awk '{print $1"\n"$4}' >SRR1039510_1.fa

#将@替换为>
sed -i 's/@/>/g' SRR1039510_1.fa

image.png

2.使用fa文件做NT比对，得到测序样本的NT比对信息

应用：随机收取2000reads做NT比对，可以发现测序样本是否含有其他物种污染等。

3.从gff或者gft文件中获取基因的ID与symbol对应关系，以及biotype类型

应用：ID与symbol转换本地化，不依赖于第三方工具和软件包，并可以根据biotype类型区分mRNA，lncRNA以及miRNA等信息。

hisat2比对

image.png
image.png

# 进入参考基因组目录
cd ~/database/genome/Ensembl/Homo_sapiens/GRCh38_release95

# Hisat2构建索引
hisat2-build Homo_sapiens.GRCh38_release95.genome.fa Homo_sapiens.GRCh38_release95.genome

# 输入输出定义文件夹
index=~/database/genome/Ensembl/Homo_sapiens/GRCh38_release95/Homo_sapiens.GRCh38_release95.genome.
inputdir=~/project/Human_16-Asthma-Trans/data/cleandata/trim_galore
outdir=~/project/Human_16-Asthma-Trans/Mapping/Hisat2

image.png

# 单个样本比对
hisat2 -p 3  -x  ${index} -1 ${inputdir}/SRR1039510_1_val_1.fq.gz -2 ${inputdir}/SRR1039510_2_val_2.fq.gz -S ${outdir}/SRR1039510.Hisat_aln.sam

image.png 第一行：fastq有多少reads，双端的reads：正端和反端各测了一遍 image.png
image.png

# sam转bam
samtools sort -@ 3 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

# 对bam建索引
samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai

# 多个样本批量进行比对，排序，建索引
# Hisat.sh内容
cat /teach/project/Human-16-Asthma-Trans/data/rawdata/sra/sampleId.txt | while read id
do
    echo "hisat2 -p 3 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log  | samtools sort -@ 5 -o ${outdir}/${id}.Hisat_aln.sorted.bam -  && samtools index ${outdir}/${id}.Hisat_aln.sorted.bam ${outdir}/${id}.Hisat_aln.sorted.bam.bai"
done >Hisat.sh

# 提交后台运行
nohup sh Hisat.sh >Hisat.log &

# 统计比对情况
multiqc -o ./ SRR*log

subjunc比对

image.png

# 进入参考基因组目录
cd /teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95

# subjunc构建索引
subread-buildindex -o Homo_sapiens.GRCh38_release95.genome Homo_sapiens.GRCh38_release95.genome.fa

# 输入输出定义文件夹
index=/teach/database/genome/Ensembl/Homo_sapiens/GRCh38_release95/Homo_sapiens.GRCh38_release95.genome
inputdir=/trainee2/Oct12/project/project_test/data/cleandata/trim_galore
outdir=/trainee2/Oct12/project/project_test/Mapping/subjunc

# subjunc单样本比对
subjunc -T 3 -i ${index} -r ${inputdir}/SRR1039510_1_val_1.fq.gz -R ${inputdir}/SRR1039510_2_val_2.fq.gz -o ${outdir}/SRR1039510.Subjunc.bam 1>${outdir}/SRR1039510.Subjunc.log 2>&1

# subjunc多样本比对
cat /teach/data/airway/sra/sampleId.txt | while read id
do
    echo "subjunc -T 1 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R             ${inputdir}/${id}_2_val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 -o ${outdir}/${id}.Subjunc.sorted.bam  ${outdir}/${id}.Subjunc.bam   && samtools index ${outdir}/${id}.Subjunc.sorted.bam ${outdir}/${id}.Subjunc.sorted.bam.bai "
done >subjunc.sh

# 运行
nohup sh subjunc.sh >subjunc.log &

SAM/BAM格式

image.png image.png
image.png

• 统计比对结果

  
  image.png

# 进入比对文件夹
cd ~/project/Human-16-Asthma-Trans/Mapping/Hisat

# 单个样本
samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam

ls *.sorted.bam | while read id
do
    echo "samtools flagstat -@ 1 ${id} > ${id/bam/flagstat} "
done >flagstat.sh

# 运行
nohup sh flagstat.sh >flagstat.log &

# 质控
multiqc -o ./  *.flagstat