三个常用参考基因组数据库：
- Ensembl：www.ensembl.org
- NCBI：https://www.ncbi.nlm.nih.gov/projects/genome/guide/human/index.shtml
- UCSC：http://www.genome.ucsc.edu/
-
Ensembl官网
image
image
image
image
image
-
参考基因组注释文件
image
image
image
image
习题
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
2.使用fa文件做NT比对，得到测序样本的NT比对信息
应用：随机收取2000reads做NT比对，可以发现测序样本是否含有其他物种污染等。
3.从gff或者gft文件中获取基因的ID与symbol对应关系，以及biotype类型
应用：ID与symbol转换本地化，不依赖于第三方工具和软件包，并可以根据biotype类型区分mRNA，lncRNA以及miRNA等信息。
hisat2比对
image
image
# 进入参考基因组目录
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
# 单个样本比对
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
image
image
image
# 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
# 进入参考基因组目录
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
image
image
image
-
统计比对结果
image
# 进入比对文件夹
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
网友评论