下载Sra数据
wget https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.2/sratoolkit.2.9.2-centos_linux64.tar.gz
tar xzvf sratoolkit.2.9.2-centos_linux64.tar.gz
vi ~/.bashrc
##添加相应路径
source ~/.bashrc
下载完成后，会在你的工作主目录下生成一个ncbi的文件夹。
cd ncbi/public/sra
for i in `seq 54 59`; do prefetch SRR70278${i}; done
图片10.png
首先进入 ncbi/public/sra$ 路径下
进行解压
fasterq-dump --split-3 SRR7027855.sra -t ./mnt -e 24
for i in `seq 54 59`
do
fastq-dump --split-3 SRR70278${i}.sra -t ./mnt -e 24
done
我这里批量运行好像有问题，因此单个解压了。
下载mt-tRNA标准序列
参考//www.greatytc.com/p/ef20c4932b3b
$ wget ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/ncrna/Homo_sapiens.GRCh38.ncrna.fa.gz
$ gunzip Homo_sapiens.GRCh38.ncrna.fa.gz
$ head -100 Homo_sapiens.GRCh38.ncrna.fa #查看fa文件前100行
查看fa文件小RNA类型
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/^>/p' | awk '{print $5}' | uniq
ENST00000459236.1 ncrna chromosome:GRCh38:7:128443449:128443510:-1 gene:ENSG00000238590.1 gene_biotype:snRNA transcript_biotype:snRNA gene_symbol:RNU7-54P description:RNA, U7 small nuclear 54 pseudogene [Source:HGNC Symbol;Acc:HGNC:34150]
TAGTGTTACAGCGCTTTGAGAATTTGTCTAGCAGGCCTTTGGGTCTTTACTGGAAAACCC
CT
sed 找到所有以>开头的行
awk '{print $5}' 输出每个项目第五个的内容
uniq 命令删除重复的行
ensemble 文件中的RNA类型
图片10.png
将其中某种类型的RNA（Mt_tRNA）提出来放入单个文件中
mkdir rfam
cd rfam
touch Mt_tRNA.fa ##新建文件
cd 在返回上一级
cat Homo_sapiens.GRCh38.ncrna.fa | sed -n '/Mt_tRNA/{p;n;p}' >> ./rfam/Mt_tRNA.fa
安装 blast+
$ wget ftp://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/2.11.0/ncbi-blast-2.11.0+-x64-linux.tar.gz
##这个安装除了问题
使用$ sudo apt install ncbi-blast+ 成功安装
makeblastdb -in ./rfam/snrna.fa -dbtype nucl -title rfam.hg.snrna -out rfam.hg.hg.trna
图片9.png
计算测序文件中不同长度的reads 总数目
nohup ls ./ncbi_database/SRR7027854.sra.fastq | while read id; do less ${id} | awk '{ if(NR%4==2) print length($0) }' | sort -n | uniq -c | awk '{ OFS="\t";$1=$1;print $2"\t"$1 }' > len.txt; done &
得到len.txt 文件，如下图
图片8.png
批量文件运行，命令行如下
nohup ls out.SRR3269*.fastq.gz | while read id; do less ${id} | awk '{ if(NR%4==2) print length($0) }'
| sort -n | uniq -c
| awk '{ OFS="\t";$1=$1;print $2"\t"$1 }' > len.${id%.*}.txt; done &
建立索引
$ bowtie2-build human_tRNA.fasta mt_tRNA
$ bowtie2 -p 10 -x mt_tRNA ./ncbi_database/SRR7027854.sra.fastq |samtools sort -O bam -@ 10 -o -> SRR7027854.bam
$ samtools index SRR7027854.bam
$ featureCounts -p -T 6 -t ncRNA_gene -g gene_id -a Homo_sapiens.GRCh38.99.chromosome.MT.gff3 -o ./untsle.txt SRR7027854.bam
其中，需要根据选择ncRNA_gene，或tRNA 或exon，然后再根据所带的属性选择gene_id. 否则会报错
图片8.png
最终结果如下，22个tRNA 匹配的长度
图片9.png
获取fastq 文件某一部分的数据
$ head -20 ./ncbi_database/SRR7027854.sra.fastq
数据类型为：
图片2.png
提取所有length=23 的序列
$ cat ./ncbi_database/SRR7027854.sra.fastq | sed -n '/length=23/{p;n;p}'
图片1.png
也可以写入新文件
cat ./ncbi_database/SRR7027854.sra.fastq | sed -n '/length=23/{p;n;p}' >> ./rfam/4snrna.fa
